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Legal Statement 


This document presents work undertaken by Frazer-Nash Consultancy Ltd and funded under contract by the 
UK Government Department for Business, Energy and Industrial Strategy (BEIS). Any statements contained 
in this document apply to Frazer-Nash Consultancy and do not represent the views or policies of BEIS or the 
UK Government. Any copies of this document (in part or in full) may only be reproduced in accordance with 
the below licence and must be accompanied by this disclaimer. 


This document is provided for general information only. It is not intended to amount to advice or suggestions 
on which any party should, or can, rely. You must obtain professional or specialist advice before taking or 
refraining from taking any action on the basis of the content of this document. 


We make no representations and give no warranties or guarantees, whether express or implied, that the 
content of this document is accurate, complete, up to date, free from any third party encumbrances or fit for 
any particular purpose. We disclaim to the maximum extent permissible and accept no responsibility for the 
consequences of this document being relied upon by you, any other party or parties, or being used for any 
purpose, or containing any error or omission. 


Except for death or personal injury caused by our negligence or any other liability which may not be excluded 
by an applicable law, we will not be liable to any party placing any form of reliance on the document for any 
loss or damage, whether in contract, tort (including negligence) breach of statutory duty, or otherwise, even if 
foreseeable, arising under or in connection with use of or reliance on any content of this document in whole 


and redistribute the material in any medium or format, provided you give appropriate credit, provide a link to 


or in part. 


Unless stated otherwise, this material is licensed under the Creative Commons 
Attribution-NonCommercial-NoDerivatives 4.0 International License. You may copy 


the license and indicate if changes were made. If you remix, transform, or build upon the material, you may 
not distribute the modified material. You may not restrict others from doing anything the license permits. 


Preface 


Nuclear thermal hydraulics is the application of thermofluid mechanics within the nuclear indus- 
try. Thermal hydraulic analysis is an important tool in addressing the global challenge to reduce 
the cost of advanced nuclear technologies. An improved predictive capability and understanding 
supports the development, optimisation and safety substantiation of nuclear power plants. 


This document is part of Nuclear Heat Transfer and Passive Cooling: Technical Volumes and Case 
Studies, a set of six technical volumes and four case studies providing information and guidance 
on aspects of nuclear thermal hydraulic analysis. This document set has been delivered by Frazer- 
Nash Consultancy, with support from a number of academic and industrial partners, as part of 
the UK Government Nuclear Innovation Programme: Advanced Reactor Design, funded by the 
Department for Business, Energy and Industrial Strategy (BEIS). 


Each technical volume outlines the technical challenges, latest analysis methods and future direc- 
tion for a specific area of nuclear thermal hydraulics. The case studies illustrate the use of a subset 
of these methods in representative nuclear industry examples. The document set is designed for 
technical users with some prior knowledge of thermofluid mechanics, who wish to know more about 
nuclear thermal hydraulics. 


The work promotes a consistent methodology for thermal hydraulic analysis of single-phase heat 
transfer and passive cooling, to inform the link between academic research and end-user needs, 
and to provide a high-quality, peer-reviewed document set suitable for use across the nuclear 
industry. 


The document set is not intended to be exhaustive or provide a set of standard engineering ‘guide- 
lines’ and it is strongly recommended that nuclear thermal hydraulic analyses are undertaken by 
Suitably Qualified and Experienced Personnel. 


The first edition of this document set has been authored by Frazer-Nash Consultancy, with the 
support of the individuals and organisations noted in each. Please acknowledge these documents 
in any work where they are used: 


Frazer-Nash Consultancy (2021) Nuclear Heat Transfer and Passive Cooling, 
Study B: Fuel Assembly CFD and UQ for a Molten Salt Reactor. 
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Introduction 


This case study is concerned with modelling a fuel assembly for a Molten Salt Reactor (MSR), 
including heat generating fluids and radiative heat transfer. To provide realism and improve the 
utility of this case study, it has been developed in collaboration with Moltex Energy’. 


As for all the case studies in this series, this analysis provides a ‘worked example’ of a specific 
modelling task to illustrate the modelling approaches described in the technical volumes. Therefore, 
it is recommended that this case study is read in conjunction with the following technical volumes: 


Volume 2: Convection, Radiation and Conjugate Heat Transfer 
Volume 3: Natural Convection and Passive Cooling 

Volume 4: Confidence and Uncertainty 

Volume 6: Molten Salt Thermal Hydraulics 


Many publicly available examples of thermal hydraulic analyses are performed as part of bench- 
mark or validation studies, as demonstrated in Study A (Liquid Metal CFD Modelling of the TALL-3D 
Test Facility) and Study D (System Code and CFD Analysis for a Light Water Small Modular Reac- 
tor). In these cases the geometry and boundary conditions are usually simple and well defined, and 
comparison or validation data (in the form of both high-fidelity simulation results and experimental 
measurements) is available. However, for many industrial thermal hydraulic modelling engineers, 
this does not generally represent the situations in which they often need to work. 


Industrial simulations are built using the best information and methods available at the time, while 
recognising that there are unknowns and uncertainties that remain, and that additional confidence 
building measures would be necessary to make a high-significance decision using it. Simulations 
of this type can be used to inform design decisions, and will form part of the process of deciding 
what additional experiments or measures are necessary to build this confidence. This case study 
and Study C (Reactor Scale CFD for Decay Heat Removal in a Lead-cooled Fast Reactor) are 
examples of this type of analysis. 


The objectives of this case study focus on two areas: 


« Using Computational Fluid Dynamics (CFD) to create heat transfer and pressure drop cor- 
relations applicable to a porous medium representation of a fuel assembly. Correlations suit- 
able for steady-state normal operation, and for quasi-steady-state buoyancy driven natural 
circulation will be developed, including an assessment of their uncertainty. 


* Quantification of the uncertainty in the CFD model predictions due to uncertainty in material 
properties. There is generally a poor state of knowledge of the thermophysical properties 


1 www.moltexenergy.com 
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of molten salts — they are not known with high accuracy, and change with temperature and 
salt composition. An assessment of the extent to which variation in them affects the figures 
of merit will be demonstrated, which can be used to target where efforts to develop bet- 
ter knowledge are best placed, and to determine which properties do not merit substantial 
additional investigation. 


While the methods demonstrated by this case study use a specific MSR design as an example, 
fuel assemblies for many reactor designs could be analysed in a similar manner. The layout of 
the described reactor is broadly in line with the Moltex Stable Salt Reactor — Wasteburner (SSR- 
W) design. However, to enable open publication, various aspects of the geometry, component 
performance and analysis scenario are fictitious. Therefore, no quantitative data presented in this 
study (either input data or analysis results) should be taken as illustrating the actual design or 
performance of the SSR-W. 
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Planning the Analysis 


Reactor and Conditions to be Analysed 


The SSR-W design uses a molten salt with dissolved fissile material contained within closed verti- 
cal fuel pins, a hexagonal array of which are arranged into a fuel assembly that has a layout similar 
to that used in liquid metal cooled reactors. The assemblies are immersed in a coolant salt in a 
pool-type reactor configuration. During normal operation, the coolant salt is driven through the fuel 
assemblies and primary heat exchangers by pumps. 


Pump and primary 


heat exchanger Inner Outer 


Wall Wall 


Upper 


plenum EHRS air 


chimney 


Active 
Core 


Lower 
plenum 


Normal operation Emergency shutdown cooling 


Figure 2.1: Schematic layout of SSR—W showing expected flow directions. 


Under the conditions of a postulated Station Blackout (SBO) event (loss of off-site and on-site 
power), the reactor will ‘trip’ (fission will stop), but the loss of power will cause the forced coolant 
flow and heat removal by the primary heat exchangers to stop. The reactor will then undergo a 
transient to establish decay heat removal by natural circulation, rejecting decay heat to a passive 
Emergency Heat Removal System (EHRS) air chimney. 


A feature of the SSR—W concept design assessed here is that during normal operation, the flow of 
coolant is from the upper plenum to the lower plenum. To establish natural circulation, this direction 
must reverse. A schematic of the reactor and flows in these two conditions is shown in Figure 2.1. 
In the passive cooling conditions, the coolant circulates through the fuel assemblies in the active 
core to the upper plenum and returns to the lower plenum via bypass paths. The decay heat is 
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transferred through the inner wall to a circulating cavity, and then through the outer wall to the air 
chimney. 


Under all conditions, the available margin between the maximum fuel salt and pin material tem- 
peratures and their limits should be known, and it is necessary to demonstrate that the EHRS is 
able to meet its required performance, while maintaining sufficient margins. The flow rates and 
temperatures in the reactor under the required range of normal operation, fault and decay heat 
removal conditions for design and licensing purposes can only be practically predicted by a low- 
fidelity representation: a system code model or porous medium CFD representation of the core’. 
These techniques require correlations to determine the flow losses and heat transfer characteris- 
tics of coolant flow in the fuel assemblies. These could be determined from representative and/or 
prototypical experiments, literature correlations or CFD analysis, potentially using a combination of 
them all to provide comparison. 


At the concept design stage of a novel reactor, it is likely that sufficiently specific experimental data 
will not be available for the geometry, fluids or flow and heat transfer regime expected. Similarly, 
literature correlations, which are typically derived from experimental data, have limitations in that 
they represent idealised geometries and conditions, and would need to be combined with each 
other in ways that are not part of their derivations. CFD does not have those limitations, and can 
also use the correlations as a way to check the accuracy and geometrical dependence of the 
results, and guide what physics and flow regime to expect. Also, at the concept design stage, the 
geometry, materials and conditions expected are subject to change, so modelling by CFD provides 
a more rapid and economical tool for design development and assessment than physical testing. 


For the SBO transient, the correlations are required to be sufficiently accurate and applicable during 
three reactor conditions: 


1. Forced convection when the pumps are active. 
2. Mixed convection after natural circulation passive cooling is established. 


3. A combination of forced, natural and mixed convection in the cross-over period following 
the reactor trip as the coolant pumps coast-down to a stop and the flow direction reverses. 
The flow in the fuel assemblies will undergo a complex transient, during which there may be 
non-uniform or counter-flow in the coolant spaces of the fuel assembly cross-section. 


This case study will assess the first two of these (the third is possible as an extension of the 
techniques described, but involves significant additional modelling challenges). 


The objective of this case study is to demonstrate how to start from a position of knowledge of 
the geometry, materials, and initial expected reactor conditions that is relevant to the concept de- 
sign stage, and use CFD to create correlations suitable for a porous medium representation. This 
includes giving some consideration to their accuracy, and to the effect of uncertainty in material 
properties. 


The long duration of the transients involved, and the need to include the whole primary and secondary secondary circuits 
means that a high-fidelity CFD model is prohibitively expensive. 
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Planning the Analysis 


As part of the ongoing reactor design process, the SBO event has been simulated using a system 
code model of the reactor, where 1/3 of the reactor (one of the three primary heat exchanger 
circuits) has been modelled with the core represented by four parallel channels (Figure 2.2), rep- 
resenting four regions from the inner (CH1) with the highest power assemblies to the outer region 
(CH4) with the least. The flow losses have been chosen such that the flow distribution among 
channels closely matches the power distribution, giving a uniform exit temperature to the lower 
plenum. 
1 of 3 250 MW primary 
heat exchangers 
2380 kg/s 
550°C 650°C 


6 FAs 54FAs 96 FAs 288 FAs 
0.0254 0.2047 0.2928 0.4771 


Figure 2.2: Core split into four concentric regions of channels, showing the number of 
Fuel Assemblies (FAs) in each for the whole core, as well as the split of core power 
and flow. Nominal core inlet and outlet temperatures are 550°C (823.15 K) and 650°C 
(923.15 K) respectively. 


The SBO transient is initiated at time t = 0, where the primary pumps are ‘tripped’ due to loss of 
power; they coast-down to stop over 100s. The nuclear heating from fission is terminated imme- 
diately (by the insertion of control rods) and a reducing decay heat is applied. This transient is 
shown for Region 2 (CH2) in Figure 2.3. The flow reverses direction at approximately 110s, and 
the upper plenum temperature exceeds the lower plenum temperature at approximately 130s. This 
case study only considers the core after this point, where natural circulation conditions have been 
established — three instants in the transient were chosen to analyse: 150s, 250s and 500s. 


Using these system code predictions as inputs, simulations can be performed of a fuel assembly 
for two types of reactor condition: 


1. Steady-state calculations for normal operation for the four core regions, applying constant 
boundary conditions. 
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Figure 2.3: Flow in CH2, core power and upper and lower plenum temperatures during 
the SBO transient. Chosen times for simulation marked. 


2. Quasi steady-state calculations of natural circulation in CH2 at the three chosen times, ap- 
plying boundary conditions that are assumed to not change significantly over the timescales 
of the flow (hence ‘quasi’ steady). 


These seven conditions are intended to capture the response of the fuel assemblies to variations 
in power and flow. The performance of the molten fuel salt is also expected to be sensitive to 
variations in: 
¢ Thermophysical and radiative heat transfer (optical) properties of the fuel and coolant salt 
and of the cladding. 
* Fuel and core geometry. 
The effect of the first of these will be explored by applying Sensitivity Analysis (SA) and Uncer- 


tainty Quantification (UQ) techniques. However, these techniques could be readily extended to 
encompass geometry variations and applied to system code simulations too. 


Moltex Fuel Assembly 


The fuel assembly design under consideration is shown in Figure 2.4 with the expected natural 
convection circulating flow pattern of molten fuel salt inside the fuel pin shown in Figure 2.5. 


The fuel pins are arranged in a triangular array within a hexagonal fuel assembly. Kakag et al. 
(1987, Chapter 7) provides a description of the geometry, frictional losses and heat transfer in 
this kind of array. The main dimensions are the pitch of the spacing of the pins, p, the pin outer 
diameter, d and the cladding thickness, t, as shown in Figure 2.5. 


The flow in the spaces between pins in the array is characterised largely by, p/d, the pitch to diam- 
eter ratio. The hydraulic diameter, D;, of the coolant flow path between the pin array is calculated 
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am 
‘ 


Top of assembly, pins have an upper 
gas space above the molten fuel salt 
and a vent for gaseous fission products 
(gas space and vent not modelled) 


Single fuel pin with wire wrap 
(wrap not modelled) 


Active section of fuel Base of assembly with representative lower 
assembly with hexagonal nozzle and fuel pin support plate with flow holes 
array of fuel pins and wrapper (nozzle and support not modelled) 


Figure 2.4: CAD model of the fuel assembly concept design, used as a source of dimen- 
sions and to visualise included and omitted features. 
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Volumetric Heat Source 


Fuel Pin (Solid) 


Coolant Salt (Fluid) 


Fuel Salt (Fluid) 


Figure 2.5: Fuel assembly pin array components and dimensions. 


Array dimensions (m) 
p 0.012 
d 0.01 
t 0.0003 
h 0.010392 
L 1.6 

D;, 0.0058783 
hquct 0.10046 


Array properties 


p/d 1.2 
Nring 9 
Npin 271 
Duct and gap dimensions (m) 
tduct 0.001 
Wgap 0.0031 


Optical path lengths (m) 

2h—d 0.010785 
p—d_ 0.002 

diye /2 0.0047 


Table 2.1: Fuel assembly geometric properties. 
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from the subchannel bounded by the triangle in Figure 2.5, where the coolant flow area is given by 


V3 1 » 
A=-—p* ——d 
4P 8 
The height, h, of the subchannel is 
es 


and the centre of the subchannel is at (2/3)h from the fuel pin centre (marked with a cross). The 
fuel assembly has 9 rings of pins around the centre pin, giving a total of 


Npin = 3Nring (Mring + 1) +1=271 


pins in each fuel assembly cross-section. The inner wall of the duct? is placed at the centre of the 
outer subchannel, so is at a distance of 


duct = (Nring + 2/3)h 


from the array centre. The total cross-section of the assembly is bounded by the duct: 


6 
Aassembly = pa hues” 


The flow area available to coolant flow is given by 


TT 


ing 


A coolant = Aassembly a Np 
The active length (containing fuel salt) of the fuel pins is L = 1.6m. The dimensions of the array 
are given in Table 2.1, along with other derived dimensions that are required later. 


Results of Importance: The primary considerations for analysing the fuel assembly are to pre- 
dict the temperature of the fuel salt and cladding material under all credible flow and heat flux 
conditions. 


* The cladding is a barrier preventing the fuel salt from mixing with the coolant which would 
contaminate the reactor, and so ensuring its integrity is a primary consideration. Damage and 
failure mechanisms can be non-linearly sensitive at elevated temperatures, and so knowing 
the maximum temperature of the cladding is typically required. 


« The duct also requires its structural integrity to be assured, and could be susceptible to ther- 
mal stresses, so the distribution of heat transfer and temperature around its circumference 
and along its length will allow this to be assessed. 


« Having a negative reactivity temperature feedback coefficient is an important operational and 
safety feature for MSRs, so knowledge of the temperature of the fuel salt is important from a 
criticality perspective. 


¢ The fuel temperature also affects the retention of gaseous fission products, and, in the most 
extreme cases, boiling of the salt must be avoided. 


2 The words ‘duct’ and ‘wrapper’ will be used interchangeably. The wrapper surrounds the whole fuel assembly, and is not to 
be confused with a ‘wire wrap’ around each pin. 
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The range of credible conditions depends on the resistance to flow in the assemblies (as well as 
the driving forces and losses in the rest of the core) and the local heat transfer between the coolant 
and the cladding. Therefore, characterising these conditions provides a direct route to assessing 
the fuel and cladding temperatures that have safety implications. 


Extent and Boundary Conditions for Fuel Assembly CFD Model 


To develop a CFD model of the fuel assembly, it is necessary to consider how the outputs from the 
simulations will be used, and what is the appropriate level of simplification that should be adopted 
in the geometry and boundary conditions. Here the target application for the results is a reactor 
scale simulation with the fuel assemblies represented by porous medium, but the results could also 
be easily converted into inputs for system code models. 


Several choices for modelling an overall core layout are described in Section 4.1 of Study C (Re- 
actor Scale CFD for Decay Heat Removal in a Lead-cooled Fast Reactor). The approach of resolv- 
ing each assembly individually, with their wrappers represented as zero thickness walls, and the 
inter-wrapper gaps also represented by porous media, has been adopted in Study C to give bet- 
ter fidelity for the prediction of core flows. The porous MSR fuel assembly model developed here 
demonstrates a compatible approach. The use of per-assembly thin walls means that the porous 
models of the fuel and inter-wrapper gap do not need to apply any specific modelling to inhibit 
flow in the across-core direction, which would otherwise be needed to account for the effect of the 
ducts. 


The inclusion of both fuel and coolant requires a Conjugate Heat Transfer (CHT) approach to be 
adopted for the fuel pin and duct, allowing the temperature distribution to be simultaneously com- 
puted in the coolant salt, fuel pin and fuel salt. Beyond this, there are a number of considerations 
involved in choosing what to model, some of which require their implications to be clearly an- 
ticipated before starting the modelling process. For others, the impact will only emerge from the 
simulation results, or exploratory precursor simulations. 


* Should the model include inlet and outlet nozzles of the assembly? 
* Does the full active length need to be modelled or would a shorter section be sufficient? 
* Should details of the fuel pin at each end be included? 


* Is modelling a whole assembly cross-section necessary, would a 1/6" sector of the hexago- 
nal fuel assembly, or even a single subchannel be sufficient? 


* Would a 2D slice be useful in determining transfer properties? 


¢ How to account for wire wraps or spacer grids? 


The design of a core and fuel to produce a desired flow distribution involves many considerations 
that combine the effects of the design of inlet and outlet nozzles to the fuel and above and be- 
low core structures. The flow openings and restrictions at the side and base of the core can, for 
example, be intended to control inter-wrapper flow and provide routes for flow under decay heat 
removal conditions. Because of this, a CFD model of a single fuel assembly cannot, on its own, be 
expected to be predictive of the core flow rates based on above and below core pressures, unless 
all these additional features are accounted for (these features can be represented in a system code 
or reactor scale CFD model, however). Therefore, to focus this work on the characterisation of the 
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fuel assembly, the flow rates for both forced and natural circulation conditions will be prescribed 
for the CFD models considered here. As long as the flow rates, temperatures and heat fluxes are 
sufficiently representative of the range of core conditions that will be encountered, the correlations 
derived will be applicable to the fuel assembly flows encountered in a reactor scale model. 


A significant simplification in this case study is the omission of the wire wraps from the fuel and 
the treatment of the pins as a bare triangular rod array. This simplifies the meshing and analysis 
of the pin array for demonstration purposes within the case study. The inclusion of wire wraps in a 
more realistic case would make some of the exploitation of simple symmetry or single subchannel 
models shown here less straightforward or not possible. Volume 5 (Section 2.3.1) discusses the 
modelling of wire wrapped fuel in the context of Liquid Metal-cooled Fast Reactors (LMFRs). Fuel 
pins without wire wraps typically have spacer grids, which can also be assessed with CFD (CSNI, 
2013). 


Omitting the wire wraps impacts on the losses and heat transfer to the coolant, but should not 
significantly affect the fuel salt modelling inside the fuel pins. The inlet and outlet nozzle structures 
will also be omitted, focusing this work on the performance of the coolant and fuel salt flow and heat 
transfer in the active length. When using the correlations developed in this study for a reactor scale 
simulation in a system code or porous medium CFD model, the effect of the inlet and outlet regions 
would need to be added on as additional losses. Separate or combined CFD model studies, similar 
to those described in this study could also be used to characterise those components. 


It would only be necessary to resolve the whole cross-section of an assembly if any asymme- 
try was to be represented, potentially arising from partial blockage or non uniform conditions for 
inter-wrapper flow, temperature or heat load in adjacent assemblies. Otherwise, a 1/6'" (or 1/12") 
symmetric sector can be used. As a further simplification, it is expected (but not guaranteed) that, 
except in the near-duct region, the behaviour of each fuel pin and coolant subchannel will be sim- 
ilar, so a model of a single coolant subchannel should be able to determine the axial behaviour. 
The use of symmetry planes for the fuel pins in a single subchannel model will force a symmetric 
recirculation pattern within them, which may or may not be physical. A model with a whole 1/6" 
sector does not force this constraint on most of the fuel pins, and so if this behaviour is present (or 
more precisely, if the CFD model can predict it), it would be visible. 


A 1/6'" sector model also gives the ability to assess the effect of the combination of thermal ra- 
diation and conduction on the effective value of thermal conductivity in the transverse direction. 
However, it is also possible to do this with a 2D slice. 


For a single subchannel or 1/6'" sector model, it would be possible to use a part-length model 
(potentially periodic in that direction) instead of the full length of the fuel. If the fuel in the pins was 
solid and static, then that would be more likely to be a practical option, but because the fuel salt 
circulates and is bounded by a solid wall at the pin base and a free surface at the top, this process 
would be expected to be subject to significant distortion if not resolved over its full height. In fact, 
the recirculating behaviour is expected to be sensitive to the heat transfer at either end, and so an 
entrance length section will be included in the model in the coolant passages to move the boundary 
away from the top and bottom of the fuel salt. 


The conclusion from these considerations is that it is worth creating a 1/6" sector of the full active 
length of the fuel (plus some entry region), resolving the full 3D physics. The results from this will 
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provide insight into what approximations are appropriate, and will serve as a reference to compare 
a 1/6" sector 3D porous model to. The porous model properties will be derived from a single 
subchannel full length model, with additional information on transverse properties derived from a 
2D slice. 


Inter-wrapper Gap: The fuel assembly duct is 1mm thick, and to complete the definition of a 
1/6'" sector of a fuel assembly, an inter-wrapper gap (distance between outer flat faces of the fuel 
assembly ducts) needs to be defined. The system code results do not include inter-wrapper flows, 
and the core flows in the concept design are not finalised, so at this stage in the analysis cycle, 
the decision was taken to match the flow characteristics in the inter-wrapper gap to the fuel bundle 
(this can be updated in future analysis as the design evolves). One of the drivers for this design 
choice is minimising thermal stresses on the duct, and so the decision to promote equal flow on 
the inside and outside of the duct is intended to mitigate larger temperature differences across it. 


By comparing literature correlations for turbulent and laminar frictional losses in a triangular array 
and a duct with infinite parallel plates, the value of wzap = 3.1mm was determined, where the 
hydraulic diameter is 2w,., = 6.2mm, similar to the 5.88 mm hydraulic diameter of the triangular 
array. The same inlet velocity and temperature will be prescribed at the entrance to the gap as 
entering the fuel bundle. 


To provide confidence that this is a reasonable choice to make for this dimension, other hexagonal 
fuel assembly designs were compared. This design detail is not readily available for many prototype 
reactors, but it was ascertained that: 


« In the MYRRHA reactor design and accompanying inter-wrapper heat transfer test rig, the 
inter-wrapper gap is 3 mm (Uitslag-Doolaard et al., 2019). 

« In the CEFR reactor design, the inter-wrapper gap is 2mm (Wang ef a/., 2020). 

« In the PRANDTL-DHxX rig the gap is 7mm, although that may be larger to facilitate experi- 


mental instrument access (Kamide et al., 2001). 


These references support the choice of the inter-wrapper gap as being realistic in terms of clear- 
ances for fuel handling, but it must be noted that they are for LMFRs (where the fuel pin arrays are 
broadly comparable, but have slightly smaller diameter pins, that are more closely spaced), which 
have significantly different heat transfer characteristics associated with the low Prandtl number 
coolant. 


Sources of Input Data 


A range of geometry, properties and operating conditions are necessary to carry out the simulations 
planned in this case study: 
* The geometry of the fuel assembly and fuel pins. 


¢ Material properties for the fuel salt and coolant (including optical properties) and the cladding, 
as well as an indication of the expected uncertainty in them. 


* The temperature and flow conditions during forced and natural circulation conditions. 


¢ The fission power (for normal operation/forced circulation conditions) or decay heat load 
(natural circulation conditions). 
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Other than the material properties, this information is available from the CAD model of the fuel 
assembly (Figure 2.4) and system code predictions. The material properties for the specific salts 
and metals have been assembled from the open literature. This information was provided by Moltex 
Energy as representative of an SSR—W concept design. 


Flow Conditions 


The flow conditions and heat source per fuel assembly in forced and natural circulation conditions, 
along with the coolant inlet temperature (taking into account the flow direction) from the system 
code results are given in Table 2.2. From this, the inlet velocity to use as a CFD boundary condition 
and the mean volumetric heat source in the fuel salt were calculated. The inlet velocity is a mean 
value to be applied over A¢oolant at the density corresponding to the inlet temperature. 


W (kg/s) Q (kW) Inlet T (K) Inlet U (m/s) vo! (MW/m?) 

Forced circulation 

CH1 30.23 3175 824.3 0.7552 105.5 
CH2 26.57 2843.1 824.3 0.6638 94.48 
CH3 21.78 2287.5 824.3 0.5441 76.02 
CH4 11.83 1242.4 824.3 0.2955 41.29 
Natural circulation (CH2) 

150s -2.404 79.037 845.1 0.06044 2.627 
250s -1.595 72.462 842.2 0.04006 2.408 
500s -1.379 60.014 845.8 0.03468 1.994 


Table 2.2: Per assembly flow, heat load and inlet temperature, and derived inlet velocity 
and mean volumetric heat source. 


Nuclear Heat Source 


The energy input from fission or radioactive decay will be applied as a prescribed heat source. 
This would typically be supplied from a neutronics/reactor physics modelling tool that would take 
into account the fuel temperature, composition and irradiation (burnup). The energy source will 
be applied to the fuel in each operating condition as a volumetric heat source, uniform over each 
fuel pin cross-section, but able to vary along the fuel assembly length. Nuclear heating in the pin 
material could also be included, but has not been in this case. 


Under normal operation, there is an axial power profile along the fuel caused by the distribution 
of neutrons, similar to solid fuelled reactors (although modified by the presence of mobile Delayed 
Neutron Precursors (DNPs)). The profile applied in this case is shown in Figure 2.6, and is used 
as a multiplier to the mean volumetric heat source (q,./) for the forced cases in Table 2.2. 


For solid fuelled reactors, there is a similar axial profile in the decay heat in the fuel, representing 
distribution and concentration of fission products. However, because the fuel salt is molten and 
mobile in this case, the fission products become sufficiently well mixed that a spatially uniform 
decay heat source can be used soon after fission ceases. Therefore, no profile is applied under 
natural circulation conditions, and the mean value is used. 
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Figure 2.6: Axial power profile for normal operation (forced circulation). 


Material Properties 


The material properties used are defined as polynomials (or piecewise polynomials) as functions 
of temperature (in K) to be easily included in the chosen CFD code. The exponential functions 
for viscosity have fifth order polynomials fit to them. The properties are plotted over the expected 
range of temperatures in Figure 2.7. The salt k and c, properties are estimated using the methods 
described in Volume 6 (Section 2.3.2). 


The source references for experimental evaluation of salts, where available, have data that is mea- 
sured at temperatures representative of the expected reactor coolant temperatures, but there are 
not well defined indications of the temperature range of validity or applicability of them. 


Similarly, there are expected to be substantial uncertainties in the salt properties. There are un- 
certainty estimates available for these properties in other common salts (for example in Romatoski 
and Hu, 2017 and Richard et a/., 2014). These uncertainty estimates are used to perform UQ (Sec- 
tion 4.3). The cladding material is a more typical and well-studied material, so its properties have 
lower uncertainties and a better definition of the temperature validity range. 


Fuel Salt Properties: The fuel salt used is KCI-UCI, (45-55 mole % eutectic) 


Density (Desyatnik et al., 1975): 
p = 4344—1.015T (kg/m?) 


Viscosity (Desyatnik et a/., 1975, converting to exponential form for consistency and to facilitate 
comparison to other salts): 


p= 1Q74:272+1644/T — 5.3456 x 1075 @3785.4/T (kg/m s) 


Specific heat capacity, estimated using the Dulong and Petit method and the molecular weight of 
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Figure 2.7: Thermophysical properties of fuel and coolant salts and fuel cladding. 
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the salt constituents, which does not provide any guidance on temperature dependence: 
Cp = 465.8 (J/kg K) 
Thermal conductivity, estimated using Khoklov’s equation and the molecular weight of the salt: 


k = —0.1963 + 0.0005T (W/mK) 


Coolant Salt Properties: The coolant salt used is KF—ZrF, (58-42 mole % eutectic) 


Density (Richard et al/., 2014): 
p = 3658 —0.887T (kg/m?) 
Viscosity (Williams, 2006, Richard et al., 2014): 
p= 1.59 x 10-*e7479/T (kg /ms) 
Specific heat capacity (estimated using the Dulong and Petit method, Williams, 2006): 
Cp = 1051 (J/kg K) 
Thermal conductivity (estimated using Khoklov’s equation, Williams, 2006): 


k = —0.032+0.0005T (W/mK) 


Cladding Properties: The cladding and duct are made from HT9 martensitic stainless steel. 


Density (Gelles, 1996): 
p= 7800 (kg/m?) 


Specific heat capacity (Yamanouchi et a/., 1992): 


T80 +500 T <800 


(J/kg K) 
3(T-800) + 550 T > 800 


Cp = 


Thermal conductivity, (Leibowitz and Blomquist, 1988): 


17.622 + 2.42 x 10-27 — 1.69 x 10-57? T < 1030 
k= (W/mK) 
12.027 + 1.218 x 10-2T T > 1030 


HT9 has been associated with advanced reactors for a number of years, and data for more of its 
properties is available in collated forms (Hales ef a/., 2016, Hofman et al., 2019). 
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Thermal Radiation Properties: The optical properties for the coolant and fuel salts in ques- 
tion are not well known, and are expected to be subject to change with the presence of fission or 
corrosion products. To account for this uncertainty, the approach taken was to bracket the range 
of possible effects by choosing absorption coefficients, x, that provide high and low optical thick- 
nesses over representative radiation path lengths in the coolant and fuel salts, ranging from nearly 
transparent to nearly opaque, as shown in Table 2.3. The three path lengths are shown in Table 2.1, 
representing the longest in-plane distance, from one fuel pin across two subchannels to another 
pin (2h — d), the shortest distance between pins (p — d) and from the centreline of the fuel to the 
cladding inner surface (dfy./2 = (d — 2t)/2). 


The chosen values of « are consistent (the same order of magnitude) with the values shown in 
Volume 6 (Section 2.5.2). A nominal value of k = 300 m~! was chosen for the fuel and coolant, to 
give balanced contributions from emission/absorption and across the gaps transmission, and the 
effect of variation between higher and lower values will be assessed. No wavelength dependence 
was implemented for the absorption coefficient, given the uncertainty assumed in its magnitude. 
A refractive index of 1.3 was chosen for both fuel and coolant salts, with no scattering, and a 
surface emissivity of 0.8 on all radiating surfaces. These are initial representative values, that can 
be revisited once the importance of radiative heat transfer and absorption coefficient have been 
ascertained. 


k (m~+) K(2h— d) k(p — d) Kdfuet /2 
100 1.08 0.2 0.47 
300 3.24 0.6 1.41 
1000 10.8 2 4.7 
3000 32.4 6 14.1 


Table 2.3: Optical thickness of representative radiation paths for assumed absorption coefficients. 


Thermal Hydraulic Phenomena of Importance 


A number of phenomena are expected to be present in the fuel assembly that are not necessarily 
typical in other flows. 


Heat generating fluid: The main differentiation between modelling the SSR—W fuel assembly and 
one suitable for a water or liquid metal cooled reactor is the presence of the molten fuel salt; 
it will be represented as a heat generating fluid and its motion and heat transfer will also be 
predicted in the CFD model. 


Participating radiation: Salts can be semi-transparent to infrared radiation, and so including ra- 
diative heat transfer in the CFD model will be necessary, treating the salts as participating 
media, that absorb and emit radiation, as well as transmit it. 


High Prandtl number fluid: The viscosity of the salts leads to moderately high Prandtl numbers, 
particularly in the coolant at the temperatures expected. This will lead to thin thermal bound- 
ary layers in the channels, requiring sufficient cells to resolve them. 


Variation in properties: The range of temperature that could be present and the strong depen- 
dence of properties on temperature, particularly density and viscosity, mean that the effect 
of temperature variation must be considered. The expected density variations mean that the 
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Boussinesq approximation is not applicable. 


While these phenomena need to be included, the first questions to ask about modelling flow in the 
fuel assembly are: are the flows laminar or turbulent? how affected by buoyancy is the flow? what 
does that mean for the expected pressure drop and heat transfer? 


Low Reynolds Number Flow and Buoyancy Effects in Coolant 


From the flow and heat source conditions in Table 2.2, it is possible to estimate the expected mean 
coolant outlet temperature from the fuel assembly, and, by combining this flow and thermal informa- 
tion with the material properties and assembly geometry, then estimate the values of characteristic 
non-dimensional numbers to assess the relevance of these phenomena. 


Inlet Re Outlet Re Bo Ra* Gr* /Re S 
Forced circulation 
CH1 1831 2875 0.0157 50867 
CH2 1609 2546 0.0218 45549 
CH3 1319 2071 0.0347 36649 
CH4 716 1125 0.152 19906 
Natural circulation (CH2) 
150s 160.1 185.3 282 70.8 1266 
250s 104.9 128.1 394 99.2 1161 
500s 92.1 111.5 380 95.6 961 


Table 2.4: Flow regime non-dimensional numbers. Bo is only applicable to turbulent 
flows; Ra* and Gr* /Re in this context are only assessed for laminar flows. 


Kakag et al. (1987, Chapter 7) describe that there is no critical Reynolds number for laminar to 
turbulent transition in rod bundles (compared to the usually clearer transition behaviour in circular 
tubes). It provides data showing that, for a triangular array of p/d = 1.2, transition begins at 
approximately Re = 1000 and fully turbulent flow is expected for Re = 4000. The effect of buoyancy 
on transition is not considered in this data, and in a real fuel assembly, turbulence at inlets from 
nozzles, or generated by wire wraps or spacer grids will be significant. Therefore, from Table 2.4 
the expected Reynolds numbers of the coolant under forced circulation conditions in Regions 1 to 
3 (CH1, CH2 and CH3) means that it is predicted to be weakly turbulent. 


The difference between inlet and outlet Reynolds number indicates the effect of the variation in 
viscosity?. For CH4, the flow is in the transition region, which complicates its assessment. Under 
natural circulation conditions in CH2, the coolant channel flow is expected to be laminar. 


The effects of buoyancy on heat transfer in turbulent vertical flow in pipes have been studied in 
detail, and are summarised by Jackson ef a/. (1989). For buoyancy assisted flows (upwards flows 
where a fluid is heated by channel walls) there is expected to be impairment of heat transfer at 
intermediate heat fluxes. Buoyancy opposed flows (downwards flows of a heated fluid) lead to heat 
transfer enhancement. 


The Buoyancy parameter Bo can be used to assess the extent of the effect of buoyancy for vertical 


For a fixed mass flow rate in a channel, Reynolds number only depends on the channel geometry and the fluid viscosity: 
Re = WD» /Aw. 
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turbulent mixed convection flows (Keshmiri et a/., 2012) 


Gr* 


23 4 
Bo =8 x 10" sas 508 


where a heat flux form of Grashof number* is used: 


« . 9g BD* 
G — 
fe kv2 


For vertical laminar mixed convection flows®, the parameter Gr*/Re (Jackson et a/., 1989) or a 
form of Rayleigh number based on the vertical temperature gradient (analogous to uniform heat 
flux conditions, Kakag et a/. 1987, Chapter 15) 


= p’BgcpD* (dT /dz) 


R 
‘a ah 


describe the extent of the effect of buoyancy. These parameters are shown in Table 2.4, evaluated 
at the mean of the inlet and outlet conditions. 


* The values of Bo for the forced circulation conditions in CH1 to CH3, where there is buoyancy 
opposed (downward) flow indicate that there is not expected to be a significant buoyancy 
effect. The Reynolds number in CH4 are too small to use this (turbulent) correlation reliably, 
and so the higher value of Bo is harder to interpret. 


¢ For laminar, buoyancy assisted natural circulation (upward) flow conditions, comparing the 
values of Gr* /Re and Ra™ to the data in their source references indicates that only a small 
enhancement to the heat transfer is expected to be added to the values of Nusselt number 
(for fully developed flow) that would be present without buoyancy. 


The effect of buoyancy on these flows is to introduce driving forces for the flow at the fuel pin 
surface due to body forces, and as such alter the velocity profile. This affects both heat transfer, 
but also friction factor. The extent to which flow losses are influenced by buoyancy is expected to 
be small, similarly to heat transfer. 


The assessments described here use correlations and results derived for flows in vertical tubes, 
so applying them to the subchannel spaces in a vertical rod array is an approximation. 


Recirculating Buoyancy Driven Low Reynolds Number Flow in Fuel Salt 


It is not as straightforward to assess the nature of the expected self-generated recirculating flow 
behaviour in the heat generating salt in the closed fuel pin, however some guidance was found 
in the literature, dating from nuclear reactor development programmes of the 1950s and 1960s°. 
Murgatroyd and Watson (1970) describe experimental and theoretical studies of the same geom- 
etry and configuration, using water with uniform internal heat generation’ and a constant tube 


4 The parameter of Gr, /Re?-’, can also be used (Jackson et al., 1989, Kakac et al. 1987, Chapter 15), but applies a different 
form of Grashof number, Gr, = g(pp — P)D?/(pv7), based on integrated density, 6 and bulk density, pp. 

5 The comparisons are made to data for experiments with uniform heat fluxes, which is a better representation of the con- 
ditions in a fuel assembly than the alternative idealisation of a uniform temperature. It is always necessary to check which 
thermal boundary condition applies to a heat transfer scenario, and to which condition a correlation or theory applies. 

® Such as the The Los Alamos Molten Plutonium Reactor Experiment (LAMPRE, Los Alamos, 1962) 

7 Achieved by passing an electrical current through 0.4 % HCl as the working fluid. 
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outer temperature. The experimental data includes resolved axial and radial profiles of velocity and 
temperature, suitable for comparison to CFD for validation. 


The flow is non-dimensionalised by the S parameter: 


gbr Guilt r r 
S= Pr-=R. 
oe Ee 
where r, is the tube radius, and q,., is the uniform heat source per unit volume. This is equivalent 
to a Rayleigh number multiplied by the radius to length ratio. 


Watson (1971) provides further analysis and plots of the experimental data, but also shows that the 
agreement of theoretical solutions with some the data should not be relied on over a wide range 
of S. Investigations of this type of configuration and flow do not seem to have been continued in 
the open literature beyond this point, possibly because the design of nuclear reactors based on 
this configuration were not pursued. The estimated value of S for the forced and natural circulation 
conditions is shown in Table 2.4. 


¢ For the natural circulation conditions, the values of S ~ 1000 are similar to the lowest values 
used by Murgatroyd and Watson (1970) and are expected to result in stable laminar flow. 


* For the forced circulation conditions, the values of S ~ 2 x 104 to 5 x 10* coincide with 
the range where much of the experimental data is present. At these higher values of S, 
some periodic instabilities or spiral motions around the tube axis were observed, and at 
S = 1.7 x 10°, the flow appeared to be turbulent. 


From this assessment, it is expected that the fuel salt flow will be laminar, but that the full active 
length of the fuel pin is needed to correctly simulate the recirculation, and that a larger cross- 
section of assembly than a single coolant subchannel is necessary to allow any instabilities and 
secondary motions, such as a spiral to develop. It should be noted that the experimental data was 
collected using an aspect ratio of r/L = 1/50, whereas for the fuel pin (dfue:/2)/L = 1/340.4. It is 
also the case that the temperature on the pin outer surface is not uniform, and in forced circulation 
conditions, the heat generation rate is not uniform with length. 


Modelling Strategy 


Based on the considerations discussed so far, an analysis scheme as shown in Figure 2.8 has 
been devised to demonstrate a process of deriving outputs that would contribute to enabling reactor 
scale simulations with more confidence. 


This involves building parts of the solution approach and modelling separately and sequentially, 
allowing a fine-grained and incremental assessment of what works well and what each part adds. 
This is possible in the scope of this case study because the fuel assembly (considered in isolation) 
is geometrically simple with well-defined flows. For larger and more complex models (such as the 
reactor-scale model in Study C) attempting this approach from the outset would be substantially 
more costly and time consuming, with less guarantee of success, because it would be harder 
to anticipate which phenomena will matter most and will present the most difficult aspects for the 
solution. For such larger cases, the approach shown here could be applied to a number of selected 
components, once their importance was better established. 
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Figure 2.8 shows the flow of information, understanding and decision making between the various 
models: 


Precursor and auxiliary models which test the applicability of the available Reynolds-Averaged 
Navier-Stokes (RANS) models to vertical mixed convection flow (described in Section 3.1); 
provide validation and mesh resolution optimisation for an internally circulating heat gener- 
ating fluid (Section 3.2), and explore the transition to turbulence in such circulating flows 
(Section 3.4). The first two were planned in advance of the simulation starting, and the last 
one was performed part-way through the process to understand the behaviour seen. 


Sixth sector assembly resolving each fuel pin and including all of the physical processes simul- 
taneously (Section 3.3). 


Single coolant subchannel modelled over the full height of the assembly to evaluate the pres- 
sure loss and heat transfer characteristics of the bundle and the maximum fuel and cladding 
temperature (Section 3.5). 


2D slice with no flow to evaluate the effective transverse conductivity of the assembly, including 
conduction in the coolant, cladding and fuel, and thermal radiation in the salts (Section 3.6). 


Porous sixth sector assembly implementing the information derived in the subchannel and 2D 
slice models, and comparing the results of the porous model to the resolved model to assess 
the accuracy of the approximation (Section 4.1) . 


Automation of the subchannel model to build surrogate models for the figures of merit that are 
not directly available from a porous model solution (Section 4.2). This technique is then 
extended to perform SA and UQ of the single subchannel to assess the sensitivity of chosen 
figures of merit to uncertainties in the salt material properties (Section 4.3). 


The reasons why the models form this particular sequence depends partly on the learning gained 
and choices made from each, and is elaborated in detail in each section. It should be noted that 
they are not simply a sequential increase in complexity — the most complex model (resolved 1/6" 
assembly) is built relatively early in the process. 


Turbulent flow is expected in the coolant, and a RANS turbulence modelling approach has been 
selected because it is a pragmatic and cost effective option, but some inaccuracies are expected. 
To provide a significantly more accurate answer, a resolved turbulence approach such as Large 
Eddy Simulation (LES) would be required, but could be orders of magnitude more costly to solve 
per model, given the size of the domain for the resolved 1/6" assembly. Also the option to use a 
single subchannel model would not be available because symmetry boundary conditions are not 
usable, and there is no natural periodicity available. The Reynolds numbers present in CH4 mean 
that it is likely to be in laminar-turbulent transition, and so for simplicity, it will not be considered as 
part of this case study. 


When implementing the porous medium models, it is expected that a non-equilibrium approach will 
be necessary, where the solid part of the porous section exists as an object in its own right, with 
thermal mass and a separate temperature field. This means that the coolant and fuel do not need 
to be assumed or forced to locally be at the same temperature, but it also requires the heat transfer 
characteristics between the coolant and cladding to be specified. 
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and cladding temperature) 
to supplement porous or 
system code models 


Figure 2.8: The sequence and flow of insight, decision making and information between 
the various models created as part of this case study, which together produce an overall 
modelling strategy. 


Modelling Tool Selection 


The ANSYS Fluent CFD code was chosen because it is available, and familiar, to the engineer 


performing the analysis, and has the key required functionality, that is able to be easily and reliably 
combined: 


* The ability to obtain steady-state solutions in closed domains without open boundaries via 
SIMPLE based pressure-velocity solvers. 

* Easy and reliable handling of CHT solids. 

¢ Non-equilibrium porous medium models. 

* Thermal radiation modelling of participating media. 

« A text-based run configuration and solution interface that is well suited to scripting a staged 


convergence strategy and to automating solutions. 


A RANS turbulence modelling approach has been chosen, but a remaining unknown at the plan- 
ning stage is whether there are adequate turbulence models available in the selected CFD code. 
This will be assessed by a precursor model. 


The open source Dakota toolkit, developed by Sandia National Laboratories will be used to imple- 
ment the SA, UQ and surrogate model aspects. It can be interfaced to Fluent to automate runs 
without significant effort and contains all of the sampling and statistical functionality required. 
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Quality Assurance and Validation 


This study is intended to demonstrate the decision making and discovery process of starting from 
a design and choosing to applying CFD to characterise aspects of its behaviour. There are no 
sources of ‘full-effect’ validation data, but it is prudent to validate the key phenomena where pos- 
sible, and perform studies to assess whether the results are sensitive to the effects of varying 
numerical modelling options. 


In addition to this, the process involves making decisions and choices about what modelling ap- 
proaches work well and should be pursued, and which should not. In a design and assessment 
setting, these choices may persist for a prolonged period and become part of the analysis ‘cus- 
tom and practice’. This introduces two possible latent conditions that may lead to inefficient or 
erroneous future work: 


« The decisions apply to the fuel design and reactor conditions assessed at this initial stage, 
and to the capabilities of the computational tools available, all of which may change, altering 
the best options for modelling. 


* Undetected mistakes may have been made at this stage because it was exploratory work that 
was not directly nuclear safety related, and these may lead to an incorrect or unsubstantiated 
conclusion being drawn. 


Therefore, the insight and outputs derived from this work should be created on the basis of a 
rigorous and robust selection process, that can be revisited in future as necessary. To assure this, 
the geometry, mesh, CFD setup, solutions, post-processing and descriptions in this case study 
report have been subject to independent quality assurance verification. Further, the simulations 
were all performed in a scripted and automated manner, so that they can be easily re-evaluated or 
extended and compared to the existing results. This is particularly pertinent in this study because it 
includes results from many hundreds of individual CFD simulations, of a number of different types, 
which build on each other sequentially. 
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The individual elements of the modelling strategy are described in the following sections, approxi- 
mately following the sequence that they were planned and executed in (although there are several 
feedback and iteration processes that have taken part along the way). The conclusions from each 
step and the decision processes for what it meant for the subsequent stages are described in each 
step, with the modelling used to determine overall comparisons, conclusions and recommendations 
discussed in Section 4. 


The precursor and auxiliary models are not described in the same level of detail as the later fuel 
assembly models, but the details of their configuration can be found in the references discussed. 


Buoyancy Affected Vertical Channel Flow 


The presence of relatively low Reynolds number flows 
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and buoyancy forces in the coolant channels raised the 
question about whether the turbulence models available 
in Fluent are able to accurately represent these effects. 
The work of Keshmiri et a/. (2012) suggests that some 
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RANS models are able to predict the reduction in Nus- 
selt number seen at intermediate values of Bo for buoy- 
ancy assisted flow, because they capture boundary layer 


relaminarisation, and that some models do not predict Full length fuel 


pin to assess 
buoyancy 
induced 

turbulence 


this. The models shown to work are not available in Flu- 
ent, but in discussion with ANSYS, it was suggested that 
their k-w SST +¥ (intermittency) transition model may be 
able to. 


This model was not designed for internal flows, however, instead it was intended to capture tran- 
sition and relaminarisation in aerodynamic applications (Menter et al., 2015). This section shows 
the results of testing this model and also assessing the performance of more common models. 


A 0.01 m long axially periodic section of a 0.1m diameter pipe was represented as a 10° wedge 
with symmetry boundary conditions, and meshed to achieve y* < 1 in the first cell (Figure 3.1). 
The flow rate was prescribed to give a Reynolds number of 5300 and the wall heat flux is modified to 
create a range of values of Bo. The fluid used was air, applying the Boussinesq approximation with 
a fixed density of 1.205 kg/m? and constant properties giving Pr = 0.71. The details of the setup 
are given by Keshmiri (2010), along with comparison data from experiments, Direct Numerical 
Simulation (DNS), LES and successful RANS results from other CFD codes. 
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Z Velocity 
0.96 


[m/s] 


Figure 3.1: Mesh coloured by axial (z) velocity for Bo = 0.25, for buoyancy assisted flow 
using the intermittency model. 


Cases where the flow is assisted by buoyancy (‘ascending’ buoyancy in Keshmiri et a/., 2012) 
show an impairment to heat transfer at moderate values of buoyancy (approximately Bo = 0.2). 
Buoyancy opposed flow (‘descending’) always show enhancement to heat transfer from the action 
of buoyancy forces. 


A large number of cases were used to assess the performance of the models and to evaluate the 
influence of initial conditions, only a subset of which are reported here in Figure 3.2. Solutions were 
converged for values of Bo from 0.02 to 100', with values clustered around Bo = 0.2. The Nusselt 
number, Nu, and Darcy friction factor, f, for the flow was calculated for each result. An additional 
simulation was run first for each setup configuration without gravity to provide a comparison without 
the effect of buoyancy, giving the prediction in the Bo = 0 limit (forced convection). These values 
are referred to as Nug and fy, and are used to normalise the results (Nu/Nug and f/f). The 
results are also compared to the Gnielinski forced convection correlation and the Moody chart, via 
the Haaland correlation (Volume 6, Section 3.1) — these values are shown in the figure titles. The 
cases reported are: 


k-w SST assisted/opposed: the k-w SST model was used, with results for buoyancy assisted 
and buoyancy opposed conditions. 

k-w SST all options assisted/opposed: the k-w SST model with the ‘low Reynolds number cor- 
rections’, ‘full’ buoyancy terms and ¥ intermittency enabled. 

k-w SST ¥ assisted/opposed: adding only the intermittency option, using a converged k-w SST 
solution first for each case. 

k-w SST ¥ (0 initial) assisted: the intermittency option added at the start of the solution with the 
74 field set to zero (laminar) everywhere. 

k-w SST low Re correction opposed: the k-w SST model with the ‘low Reynolds number correc- 
tions’ enabled. 

realizable k - ¢ buoyancy assisted/opposed: the realizable k - € with ‘full’ buoyancy terms enabled. 


The ‘k-w SST ¥ (0 initial) assisted’ simulation produced the same result as the -y model results 


1 Only from 0.02 to 2 for some configurations, because the very high heat flux cases were of less relevance. 


25 of 103 


Performing the Analysis 


that are started from a converged k-w SST result for Bo > 0.16. However, at lower values of 
Bo, the simulations that were initialised with zero -y produce values of Nu and f significantly lower 
than Nu and fp because they erroneously predict laminar behaviour across a large portion of the 
cross-section of the channel. This indicates that the -y model is unable to correctly generate its own 
turbulence from some starting points. The same effect is seen if the solution is initialised with a 
uniform -y = 1, and is also observed when using an approach where the solution at one value of 
Bo is started from the converged result from the next higher value of Bo (the solutions being run in 
sequence without being reinitialised each time). This is interpreted as indicating that the solution 
produced by the -y model includes an element of hysteresis, or solution trajectory dependence, for 
whether or not it predicts the correct boundary region laminarisation. 


From the results of all of the cases, it can be seen that. 


All of the RANS models tested over-predict Nuo and fo compared to the correlations for forced 
convection, although the relatively low Re flow may make this comparison more challenging. 
The y model results are close to the experimental, DNS and LES results reported by Keshmiri 
(2010). 


The k-w SST ¥ intermittency model (with default parameters) is the only case plotted that 
performs well for the buoyancy assisted and opposed flows for Nu and f. For buoyancy as- 
sisted flows it predicts the relaminarisation and heat transfer impairment behaviour (reduction 
in Nu/Nug to = 0.4 at around Bo = 0.2, with corresponding reduction in f /fo) to the same 
extent as the other successful RANS and LES results (from Keshmiri, 2010) that are plotted 
for comparison with the experimental data. It does appear to be sensitive to solution strategy, 
however. It also produces good predictions of the enhancement to Nu and the behaviour of 
f for buoyancy opposed flow, compared to the available experimental and DNS results. 


The performance of the k-w SST model is in-line with the behaviour found by Keshmiri ef al. 
(2012), where heat transfer impairment was also found to be not well predicted. 


The complexity of the numerical implementation of the models, and of the physics of buoy- 
ancy generated relaminarisation behaviour makes it hard to tell (without significantly more 
detailed investigations) if the ‘hysteresis’ effect is present for true physical reasons, or is 
purely a modelling artefact. For these relatively low Reynolds numbers, delayed transition is 
known to occur in the absence of disturbances (Eckhardt et al., 2007, Mullin and Peixinho, 
2006) and has been shown to be reluctant to transition in DNS simulations (Wu et al., 2015). 


Starting from a converged k-w SST solution before enabling -y appears to be the most 
reliable and general approach for using the intermittency model in this case. 


The ‘low Reynolds number corrections’ to the k-w SST model worsen its performance at 
resolving the heat transfer impairment and increase the values of Nup. When combining this 
option with the + model, it is predominately the ‘low Reynolds number corrections’ behaviour 
that is observed. 


Adding the ‘full’ buoyancy terms to the k-w SST or realizable k - « models has little effect. 


From Table 2.4, buoyancy opposed flow in the coolant channel at Bo < 0.035 is of most relevance 
to this case study, although for flow at lower Reynolds numbers (Re < 3000). Details of the propen- 
sity for and extent of laminarisation in the rod array under these conditions is not well known, and 
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Figure 3.2: Nusselt number and friction factor ratios for buoyancy assisted (left) and 
opposed (right) cases. Results with solid lines are new simulations performed here, sym- 
bols are comparison experimental, DNS, LES or RANS results. 
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cannot be inferred reliably from these vertical tube results. The turbulent coolant flow in the fuel 
assembly that is being considered is all buoyancy opposed, which should promote turbulence, and 
so some enhancement of Nu is expected, but only by a small amount. The expected effect on f is 
not as clear, but also not substantial. 


Despite the promising performance of the + intermittency model, it will not be used in the follow- 
ing simulations. This is because, there is a risk of the solutions exhibiting the ‘hysteresis’ effect 
found, because low levels of intermittency may be transported from laminarised regions into higher 
Re sections, where they could cause problems with creating or sustaining the correct state of tur- 
bulence. This would be possible to examine and explore for a small number of simulations, but 
because these models are intended to be run hundreds of times, being sure of the correctness 
and robustness of this behaviour in each would require assessment of each, adding a potentially 
significant burden. 


Closed Vertical Tubes With a Heat Generating Fluid 


One of the aspects of the fuel assembly that is unfamil- 
iar and uncertain is the closed recirculating fuel salt with 
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internal heat generation. It is not clear whether there will 
be a stable flow within the tube and whether there will 
be one full-length recirculation, or a number of convec- 
tion cells, as found in tall parallel plate cavities (Jin and 


Validation of Sixth sector 


buoyant flow in assembly 
Chen, 1996). Bergholz (1980) suggests that there are not Bacar Be carn aeel Pecolving 
heat generating pins, cladding, 


expected to be instabilities that cause a cellular structure 


fluid fuel and duct 


to occur, and a stable upward flow at the centre of each 
pin is likely. Full length fuel 
pin to assess 
buoyancy 
induced 
turbulence 


The experiments and analysis of Murgatroyd and Watson 
(1970) and Watson (1971) provide spatially resolved val- 
idation data, insight into the expected phenomena, and 
a straightforward case to test solution methods against. 


While the design feature of molten fuel in long closed tubes is specific to the Moltex design, the 
requirement to simulate buoyancy driven recirculation of variable density fluids in a closed space, 
or regions of high heat load but low circulation, is a commonly encountered situation. The process 
shown here of finding validation data, accompanied by dimensional analysis, developing meshing 
approaches and testing the necessary resolution, and proving the numerical solution route are a 
transferable pattern and thought process. 


Another aspect that required testing in a simple case is the simulation of a closed domain con- 
taining a fluid with variable density. In reality the fuel has a gas space above it, but it would be 
unnecessarily complex to try to simulate this gas space explicitly using a multiphase approach. 
Instead, a pragmatic approximation is to apply a zero-shear wall (slip) velocity field boundary con- 
dition to the top of the domain to represent the free surface. Doing this creates a closed domain. If 
the Boussinesq approximation is used then the density of the fluid is fixed for the purposes of the 
continuity equation, and the mass in the domain is fixed. However, the density of fuel salt depends 
on temperature, and the Boussinesq approximation is not applicable to these cases with large 
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changes in temperature over the tube cross-section, and with axial distance. In this situation, the 
mass of salt in the fluid must change as the temperature solution converges, and there is no link to 
the pressure of the domain in the equation of state (for liquids which are nearly incompressible, this 
is generally a good approximation). For closed spaces containing gases, this could be accounted 
for by conserving mass and allowing the pressure in the space to solve to the physically correct 
value. 


The geometry and mesh is shown in Figure 3.3. The mesh shown contains approximately 100,000 
cells, and was optimised in wall normal and axial dimensions over four loops of iterative testing to 
give results almost indistinguishable from a 770,000 cell mesh (mainly achieved by stretching the 
axial cells to 20 mm in length and grading them towards the top and bottom). This mesh contained 
a small ‘pin-hole’ aperture/vent in the top surface that could be made either a zero-shear wall (to 
close the domain) or a pressure boundary, able to allow mass flow in either direction. 


O-grid pattern well resolved 
in the wall normal direction 
with 0.5 mm first cell height. 
Extruded through height with 
varying vertical resolution. 


Small 'pin hole' aperture 
in top at free surface to 
test solution options for 
mass conservation in a 
closed system with 
varying density. 


Graded vertical resolution from 1 mm at base 
(wall) to 20 mm over the majority of the length, 
and reduced to 2 mm at the top (free surface). 


Figure 3.3: Geometry to compare to experimental data from Murgatroyd and Watson 
(1970) and Watson (1971) (length L = 1270 mm, radius r = 25.4mm, r/L = 1/50). 


The details of the experimental external boundary conditions, heat sources or temperatures used 
are not given precisely by Murgatroyd and Watson (1970) or Watson (1971). They are either omit- 
ted or presented in non-dimensional form only. For a given reference temperature, the volumet- 
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ric heat source, gyo; required for a given value of S can be found using the expression in Sec- 
tion 2.2.4.2 by evaluating the material properties. However, the temperature to evaluate them at is 
not stated. 


The best information is that the sources refer to the Pr value for water in the range of temperatures 
used as being 3 < Pr < 8 which gives a range of temperatures of approximately 15 to 50°C. 
It is also not clear whether it is the fluid temperature or the reference (wall) temperature that is 
used to evaluate the properties, and the source of the water properties is not stated, and hence no 
assessment can be made of their accuracy. 


The experiment is designed to provide a close approximation of a constant temperature on the 
outer surface of the perspex tube containing the fluid, and so an external convection boundary 
condition was applied using a Heat Transfer Coefficient (HTC) of 100 W/m? K and an external 
temperature of 300 K. 305 K was used as a representative fluid temperature to determine the fluid 
properties and heat source. These were sourced from the IAPWS-IF972 (Wagner et a/., 2000) and 
are shown in Table 3.1. 


293.15 K 305 K 
p (kg/m?) 998.21 995.08 
k (W/mK) 0.5995 0.6177 
pw (kg/ms) 0.001002 0.0007668 
Cp (J/kg K) 4184.8 4179.5 
Pr 6.991 5.189 
B(1/K) 0.0002069 0.0003191 
voi (W/m?) for S = 8900 1790.8 950.3 
voi (W/m?) for S = 84300 16963 9001.3 


Table 3.1: Water properties and required volumetric heat sources for two values of S. 


The necessary values of q,.; for two example values of S used for the experiments are also in- 
cluded. It can be seen that reducing the reference temperature to 293.15 K results in a difference 
in the heat source terms of a factor of 1.9, driven mainly by the changes in 1/v? and 8. This sig- 
nificant sensitivity to operating temperature introduces a large uncertainty in the comparison of the 
CFD results to the experiment for a given value of S. 


The case was setup in Fluent, as the intended CFD code, but also in CFX because it uses a 
different numerical approach to the solution, and also to test the use of closed domains in CFX, 
which is relevant to Study C. Three setup options were explored: 


« The water density represented by the Boussinesq approximation, with a fixed value of density 
evaluated at 305 K, or a variable density implemented as a polynomial with a dependence 
on temperature only. 


* The small aperture was either open (able to exchange mass in either direction) or closed (a 
zero-shear wall as part of the free surface). 


« Pressure-velocity coupling using either an algorithm based on a segregated solution using 
the SIMPLE algorithm, or a coupled solution. The SIMPLE algorithm is not available in CFX. 


2 Using an open source MATLAB implementation available from github.com/mikofski/[APWS_IF97. 
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All solutions were configured to run in steady-state mode — for SIMPLE this is a true steady- 
state solution, for coupled solutions, the numerical time step is manipulated to reach steady- 
state. 


The solutions attempted are shown in Table 3.2, which were all run with laminar settings. 


« Using the Boussinesq approximation, the cases solve well in both codes, and no opening is 
needed. These also form good starting conditions for variable density cases. In Fluent when 
allowing a variable density, running with a SIMPLE based solver was successful with and 
without using an opening — the changes in system mass were absorbed by the algorithm 
when enforcing continuity. In Fluent when using variable density, a converged steady-state 
solution was not obtained using the coupled solver with the settings tested’, whether includ- 
ing an opening or not. 


« In CFX, for variable density without an opening, the solution sometimes converged, but the 
velocity field was occasionally subject to sudden large excursions. 


« In CFX, for variable density with an opening, the solution converged, but chose the wrong 
numerical timescales by default, so converged very slowly. A manual timescale estimate was 
necessary (10 to 100 times bigger) to obtain rapid convergence. The Boussinesq solution 
estimated the timescale correctly, and so was a good guide for what to apply. 


Density Aperture P-U coupling Fluent CFX 
Boussinesq Closed SIMPLE Stable 
Boussinesq Closed Coupled Stable 
Variable Closed SIMPLE Stable 
Variable Open SIMPLE Stable 
Variable Closed Coupled Unstable Unstable 
Variable Open Coupled Unstable Stable 


Table 3.2: Solution attempts for closed domain. 


Monitoring the mass of fluid in the domain (and the mass flow rates across the opening, where 
present), the minimum and maximum vertical velocities, maximum temperature, and wall heat flux 
were used as guides to assess convergence and energy conservation, in addition to assessing the 
solution residuals. 


To obtain any converged solutions with a variable density, it was found that a good estimate of 
operating density was necessary. The operating density, p., redefines the pressure field for the 
solution to be P’ = P — popgz, where z is the co-ordinate in the direction of gravity. This removes 
a reference hydrostatic contribution’? from the pressure field for numerical reasons. It improves 
convergence by separating the small pressure differences generated by buoyancy forces and flow 
from the large hydrostatic variation®. A value of pop = 995.08 kg/m? was used. 


3 Including with or without the ‘pseudo-transient’ setting. 

4 It cannot be said to remove the hydrostatic contribution unless the density is constant. When density varies, the hydrostatic 
pressure is the result of the integral of the density field above a point. 

5 This is true even when using a double precision solution. All simulations reported in this case study use double preci- 
sion, and it should be the default for all CFD solutions, except in specific situations of relatively simple physics, such as 
incompressible external aerodynamics, and where the usable mesh sizes are limited by available memory. 
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Figure 3.4: Comparisons of CFD predictions of non-dimensionalised vertical velocity 
and temperature to the available experimental data (symbols). 


The solutions produced by the two CFD codes were very similar, where they could both be ob- 
tained. The comparisons between codes and the assessments of the iteration of mesh resolution 
were made by comparing the details of the axial variation of centreline temperature and vertical 
velocity. 


This exploration shows that using Fluent with a steady-state SIMPLE solver, using an operating 
density and starting a simulation from a Boussinesq condition before switching to a variable density 
fluid provides a reliable solution, without needing an open aperture. CFD results produced using 
this approach were compared to the experimental results combined from both Murgatroyd and Wat- 
son (1970) and Watson (1971). They were compared using radial profiles of non-dimensionalised 
axial (vertical) velocity and temperature at mid-height (x = 0.5, where x is the non-dimensional 
axial distance). The non-dimensional quantities are defined as 


_ Ur? 


Tol 


and 
k(T —T,) 


vol re 


where 7, is the local inner wall temperature at that height. The results are shown in Figure 3.4, 


t= 
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where the comparisons are split into those at high and low S$ for clarity. The t values are plotted 
relative to their centreline value to make the result coincide at r = 0 to emphasise the differences 
in shape. The agreement is convincing, showing the peaking effect in the temperature field that 
coincides with the radius where zero velocity occurs. More detailed numerical comparison of the 
level of agreement has not been pursued because of the sensitivity of S for a given heat load to 
the operating temperature. 


Aside from demonstrating an initial validation comparison, giving confidence that CFD is able to 
predict these laminar flows well, performing this precursor simulation was a useful step for a num- 
ber of reasons. It allowed a solution and monitoring strategy to be developed, and practised with 
a simple and small model, which could then be applied to the resolved fuel assembly case. It also 
provided guidance on the appropriate mesh resolution for the resolved model, and gave confidence 
that an aperture does not need to be included, simplifying the meshing and setup process of the 
larger model. 


Resolved Rod Assembly 


The largest and most complex model solved in the case 
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study was the 1/6'" sector resolved geometry, containing 
all of the components and physics simultaneously. 


For brevity, only one result for forced circulation condi- 
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tions (in CH2) and one for the natural circulation condi- 
tions (at 250s in CH2) were simulated (these will be re- 
ferred to as the ‘baseline’ conditions), although the meth- 


ods described below would allow any or all of the flow 
Full length fuel 


pin to assess 

buoyancy 
induced 

turbulence 


conditions to be evaluated. 


Geometry and Mesh 


The mesh for the resolved mesh was created using ANSA. This tool was selected because it is 
powerful and flexible for producing quad dominated patterns on surfaces, and then using them to 
create extruded meshes with axial resolution control. It also has the facility to ‘link’ faces in the 
geometry, so that the cross-section for a single fuel pin, coolant subchannel and piece of side 
symmetry could be created and then replicated. Changing the mesh on any one of these linked 
faces changed it on all of them automatically, allowing modifications to the mesh to be made rapidly 
and reliably. 


The mesh was built from an analytical description of the geometry, rather than retaining any of the 
CAD model geometry itself. The CAD model was imported and overlaid to check that the major 
features were the correct size and in the correct place. The locations for the repeating units for 
linked faces are shown in Figure 3.5. These points were generated analytically and imported into 
ANSA as a text file of positions to minimise manual transcription errors. 


The mesh pattern chosen for the cross-section is shown in Figure 3.6. The fuel pins were meshed 
with an adapted O-grid pattern based on six segments (to allow the centre and symmetry plane 
pins to be identical to those in the main assembly). The mainly quad face pattern was allowed to 
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Figure 3.5: Locations of the centres of fuel pins (red) and subchannels (blue) to provide 
replication points for linked units of mesh, and the location of the duct inner and outer 
surface, and centre of the inter-wrapper gap (green). 


introduce triangle faces to reduce the circumferential number of cells with reducing radius. This ap- 
proach is also used in the cladding and boundary layer region of the coolant. Without this approach, 
and using quad faces only, an excessive number of converging radial lines would be present at the 
centre of each fuel pin. The cladding and duct are included explicitly as solid regions and so will be 
solved for using a CHT approach. 


This profile was extruded over the 1600 mm active length of the fuel. 100 mm was added at the top 
to allow an entrance region for the coolant, and 10 mm was added at the bottom, representing the 
gap between the bottom of the fuel salt and the pin support plate (Figure 2.4). The resulting mesh 
is shown in Figure 3.7, where the axial and planar resolutions were guided by the optimised mesh 
from Section 3.2. The axial resolution expands to 4mm over the first 100 mm from each end of the 
active length, and then to 20 mm over the next 100 mm, so is constant at 20 mm over the middle 
1200 mm of active length. 


Other key cell size are: 


* 0.03 mm first cell height (normal to surface) in the coolant at the duct and cladding surface. 
The cells expand in layers, and have a target size of 0.16 mm in the middle parts of the 
coolant gap between fuel pins. An estimate from the expected flow rates indicated that a first 
cell height of 0.06 mm would be necessary for yt = 1 in the coolant, so the near-wall flow is 
expected to be well resolved. 


* 0.09 mm first cell height in the fuel salt normal to surface. 
* 0.3mm circumferential cell size on fuel pin surface. 


* 4 cells were used through the cladding thickness. 


The resulting mesh contained 44 million cells. This approach to meshing results in cells that are 
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Figure 3.6: Mesh cross-sections of fuel pin, coolant subchannels and duct region for the 
1/6" resolved fuel assembly. Only showing a subset of the mesh at the corner of the 
domain. 


almost entirely hexahedra, with a small number of triangular prisms and no pyramids, tetrahedra or 
polyhedra. The resulting mesh had a worst cell quality (using Fluent’s metrics, ANSYS, 2020b) of 
0.846 equivolume skew (where values > 0.95 indicate poor quality, that would badly compromise 
local accuracy or overall solution stability) and 0.231 orthogonal quality (where numbers close to 
zero indicate poor quality, and < 0.01 would lead to solution stability difficulties). Because they are 
mainly flow-aligned hexahedra, most cells in the mesh are of significantly higher quality than these 
worst values. 


The high quality cell shapes give confidence that the solution should not exhibit divergences for 
mesh quality reasons. High aspect ratio cells are present, and not just in boundary layers. This 
is due to the stretched 20 mm axial resolution necessary to limit the cell count. Because the cells 
are stretched in the flow direction, and the axial variation in velocity and temperature is slow and 
smooth, this does not cause problems for converged, stable, unidirectional flow solutions. If the 
flow was subject to recirculation in these regions, these cells could be problematic for the solution. 


A non-conformal interface was applied between the cladding and the fuel salt over 100 mm at 
either end of the fuel pin — without it, the fine axial resolution of the fuel would need to be present 
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in the mesh for the coolant channel. Using the non-conformal interface in this way reduces cell 
count and avoids unnecessary changes in refinement in the flow direction in the coolant, whilst 
allowing the flows at the ends of the fuel salt to be sufficiently resolved. The local axial temperature 
gradients present in the fuel salt at the ends of the pins, in the sections covered by the non- 
conformal interface, are not sufficiently high that the coarser axial solid (and coolant) would be 
unable resolve them. 


Case Setup 


Materials were created with polynomial dependences on temperature as described in Section 2.2.3.3 
and applied to the fuel salt, coolant salt, cladding and duct. 


Inlet and Outlet: For the forced circulation cases, the coolant inlet is at the top of the domain; 
for the natural circulation cases the inlet is at the base. The flow was introduced with a uniform 
coolant velocity and temperature at the values given in Table 2.2. Uniform inlet conditions are an 
approximation: 


* For the forced circulation flow, the inlet is 100/5.88 = 17 hydraulic diameters from the top 
of the active section containing fuel, so a substantially well developed velocity profile is ex- 
pected to exist®. 


¢ For the natural circulation flow, the coolant must pass through the holes in the support plate 
only 10 mm below the active section, and so will not be a developed profile’. 


The inlet turbulence values were obtained from a pipe flow correlation based on the inlet Reynolds 
number and hydraulic diameter (ANSYS, 2020b). The outlet for the coolant is set to zero gauge 
pressure. For simplicity, the same inlet and outlet boundary conditions were applied to the inter- 
wrapper gap as to the coolant in the fuel pin bundle. The bundle and gap fluid domains are sepa- 
rated by the duct and do not exchange flow or communicate pressure. 


Outlet boundary conditions usually require temperature and turbulence values to be specified for 
flow introduced by any backflow. However, this case is a unidirectional channel flow in the coolant, 
so no backflow could occur physically. As a solution convergence improvement, the outlets were 
set to create wall faces to prevent backflow if, at an intermediate point during the solution process, 
the pressure field was such as to generate any reversal. 


Heat Source: The heat source given in Table 2.2 was applied to the fuel salt regions of the 
domain, with the spatial distribution shown in Figure 2.6 for the forced circulation conditions imple- 
mented as an axial profile. 


Longer entrance lengths would be needed to produce fully developed velocity profiles, but these would be longer than 
the available upstream distances in the actual fuel assembly. The cladding outer surface over the entrance length is also 
approximated as being adiabatic, but the coolant environment above the top of the fuel salt would create some heat transfer, 
so there would also be a partially developing thermal profile present in reality. 

The design of the holes in the support plate is schematic only at this stage in the reactor development, and so detailed 
consideration of flow through them is not warranted. 
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Half of an inter-wrapper gap, with 
symmetry on outer face. 


Coolant inlet/outlet to top plenum for 
bundle and inter-wapper gap. 


” 


Wrapper/duct 


Inter-wrapper gap 


Sixth sector of full assembly 
representing 1600 mm active 
length containing fuel salt 
with a 100 mm space above 
and 10 mm space below in 
the coolant region only. 


Symmetry on 60° planes 


Coolant inlet/outlet to 
bottom plenum for bundle 
and inter-wapper gap. 


Top of fuel salt Non-conformal interface between 
(free surface) cladding and fuel salt. 0.4 mm 
axial resolution at free surface. 


Fuel 

Six-sector O-grid pattern with well resolved 
boundary layers extruded through height with 
varying vertical resolution. 


Coolant 
Well resolved boundary on cladding outer surface 
and mainly quad cells extruded axially. 


Cladding 

Mainly quad faces extruded to form hex cells - 

conformal with coolant, axially non-conformal with 

fuel over 100 mm section at either end. Non-conformal interface between 
cladding and fuel salt. 0.2 mm axial 
resolution at lower wall in fuel. 


Figure 3.7: Geometry and mesh of the 1/6" resolved fuel assembly. 
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Operating Conditions: The coolant outlet pressures are relative to a reference pressure of 2 bar 
applicable in the domain coolant. The fuel pins are closed, and have their pressure referenced to 
the same value of reference pressure at a fixed location at the top of the fuel pins. However, the 
absolute pressure in the domain has no effect on the solution, because the density of the fluid 
depends only on temperature, and none of the other material properties depend on it. 


The operating density was set to be the fuel salt density at a temperature estimated to be represen- 
tative of the final solved value. For the forced circulation case, this was 3126 kg/m? representing 
1200 K, and for the natural circulation case it was 3430.5 kg/m? representing 900 K. 


Turbulence: The k-w SST RANS model was applied, and the low Reynolds number corrections 
were activated. This was done, despite the findings of Section 3.1, because a stable solution proved 
difficult to obtain without them. The choice of this model is discussed further below. 


Participating Radiation: The Discrete Ordinates (DO) radiation model was chosen to represent 
the salts as participating media, and the nominal material and surface properties described in Sec- 
tion 2.2.3.3 were applied. The radiation models typically available in CFD codes are described in 
Volume 2 (Section 3.4.3.2), and the considerations necessary for participating media are described 
in Volume 6. DO is the most suitable model because it can accommodate any optical thickness, 
and can be easily extended to represent a banded spectral (wavelength dependent) absorption 
coefficient. 


It was found to be necessary to have a minimum level of angular discretisation and pixelation (AN- 
SYS, 2020a, Murthy and Mathur, 1998) to avoid spurious results being caused in the temperature 
field at the 60° symmetry planes. It was necessary to use 3 x 3 angular discretisation per octant of 
angular space, and 5 x 5 pixelation of the control angles. 


Obtaining a Converged Solution 


Obtaining a stable solution for this case was not straightforward. It was found to be necessary 
to carefully control the numerical settings, turbulence modelling and solution initialisation, as well 
as adding some aspects of the solution incrementally. At each stage, if any of these conditions 
were not appropriate, the solution residuals would increase, the coolant outlets would attempt to 
backflow over a significant fraction (>50%) of their area, introducing wall faces, and eventually 
the solution would diverge, resulting in a floating point error. The stages of the solution strategy 
described below are summarised in Tables 3.4 and 3.5. 


Initialisation: The solution initialisation is summarised in Table 3.3. In the forced circulation case, 
it was found to be advantageous to provide an initial vertical velocity, and in both forced and natural 
circulation cases it was found to be beneficial to provide domain specific initial guesses for tem- 
peratures to the cladding and fuel salt. For the turbulence fields in the coolant in both cases, the 
inlet value estimates were applied for the k and w fields. For the fuel salt in the forced circulation 
case, the turbulence in the coolant was too high, and the same low turbulence values as the natural 
circulation case were applied. Zero gauge pressure was applied everywhere. 


38 of 103 


Performing the Analysis 


Forced circulation Natural circulation 
Coolant Fuel Coolant Fuel 
P (N/m?) 0 ) 
U (m/s) (0 0 -0.05) (0 0 0) (0 00) 
k and w As inlet As natural As inlet 
Coolant Fuel/Clad Coolant Fuel/Clad 
T (K) As inlet 1200 As inlet 890/850 


Table 3.3: Initialisation settings for forced and natural cases. 


Source and Buoyancy: It was found necessary to obtain a relatively well converged solution 
for the flow and temperature fields in the coolant, to provide a stable boundary condition to the 
cladding outer surface, before attempting to solve for the flow in the fuel salt. This is achieved by 
starting the case with gravity off, and adding the heat source to the fuel incrementally, starting at 
25% for an initial period to allow the flow to develop before stepping up to 100 % and converging 
the temperatures. While gravity is off, there is no circulation in the fuel salt, so energy can only be 
transported by conduction. To limit the rise in centreline temperature, the fuel salt conductivity was 
artificially increased by a factor of 10. Gravity was then introduced in three steps, with the number 
of iterations used between each increment chosen by monitoring the behaviour of the solution to 
observe a tendency towards convergence before each change. 


A key method for obtaining a solution for the salt flow is running with the Boussinesq approxima- 
tion in the fuel salt to establish the flow and temperature fields before allowing density to vary. 
Because it is an open channel and there was less local variation in temperature, the coolant den- 
sity was allowed to vary with temperature throughout. While running with Boussinesq, the fixed 
fuel salt density was set to be the same as the operating density for each case (3126 kg/m? or 
3430.5 kg/m?). 


This solution is harder to run because there are disconnected fluid regions that are strongly affected 
by buoyancy forces, with substantially different fluid densities (e.g. 3430.5 kg/m? in the fuel vs. 
2860 kg/m? in the coolant at 900K), but only a single operating density can be specified. The 
operating density was chosen to suit the fuel because that domain is more challenging to converge. 
Similarly, it is also only possible to run a single RANS model in both domains (although it is possible 
to solve for laminar flow in a subset of domains), which may not be ideally suited to both. 


Numerical Settings: Based on the experience of the cases described in Section 3.2, the SIM- 
PLE segregated pressure-velocity coupling scheme was chosen, with body-force weighted pres- 
sure discretisation applied, as recommended for buoyancy driven flow (ANSYS, 2020b). By ob- 
serving the solver messages generated during the early solution process, it was also found that 
applying the BCGSTAB pressure stabilisation method was beneficial. 


The solutions were started with first-order upwind discretisation for the momentum and energy 
field, and were changed to second-order upwind discretisation when the solution was well enough 
converged to allow ‘full’ buoyancy (variable density) in the fuel salt, rather than the Boussinesq ap- 
proximation. The fixed densities for the Boussinesq representation were chosen to be a good esti- 
mate of the mean density that would exist in the fuel at the temperature reached at this switchover 
point, so the jump in system mass is minimised, aiding the convergence of mass continuity. 
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The Under-Relaxation Factors (URFs) for the density, turbulent viscosity, temperature and body- 
force terms were reduced from their default values of 1.0 to 0.95 for the start of the solution and 
increased to 0.98 later. The default values for other URFs (particularly 0.3 for pressure and 0.7 for 
momentum) were retained throughout. 


The thermal radiation solution is computationally expensive and only changes as a result of the 
evolution of the temperature field (which depends on the flow), so the radiation solution was only 
turned on after the temperature and flow fields had converged. In addition, the radiation fields were 
only solved every 4 flow iterations for the forced circulation case, and every 10 iterations for the 
natural circulation case. 


At the end of the natural circulation case, the turbulence modelling was turned off, and the solution 
only solved for laminar flow. Despite the flow being laminar in its final solution, it was not found 
to be possible to run the case without a RANS model active in the early stages, to damp out the 
‘spikes’ in the velocity fields that occur, particularly with the addition of gravity. 


Iteration URFs Source Gravity Buoyancy U,T order Radiation 
start 0.95 25% 0 Boussinesq 1st off 

50 0.98 T 

100 100% 

300 -1 

325 5 

400 -9.81 

600 0.98 full ow 

982* on 
1769* 0.98 100% -9.81 full on on 


Table 3.4: Sequential steps of the solution strategy for the forced resolved rod case. 


Iteration URFs Source Gravity Buoyancy U,7T order Radiation Turbulence 
start 0.95 25% 0 Boussinesq 1* off RANS 

25 0.98 T 100% 

250 -1 

275 -5 

300 0.98 -9.81 

400 full and 

1400 on 

2824* laminar 

2875* 0.98 100% -9.81 full on on laminar 


Table 3.5: Sequential steps of the solution strategy for the natural resolved rod case. 


Monitors and Convergence: The convergence of the solutions was assessed by monitoring the 
progress of the solution residuals, and also key quantities providing insight into the solution state 
and stability, and that also provide criteria for judging when a steady-state had been achieved. 
Because of the lack of features such as recirculation and separated flow, all of the solutions in 
this case study could be converged to true steady-states. This would not be the case with more 
complex geometries, where small regions of unsteadiness could prevent it. 


The following quantities were recorded at each iteration: 


40 of 103 


3.3.4 


Performing the Analysis 


Maximum temperatures: The maximum cell temperature in the bundle coolant, fuel, cladding, duct 
and gap indicate if the solution remains bounded between acceptable values, and when the 
flow and heat transfer processes have stabilised. Jumps in the value or gradient of these 
monitors were seen as the steps of the solution strategy were applied. The application of 
thermal radiation results in a substantial reduction in maximum fuel temperature, for example. 

Fuel salt velocites: The minimum and maximum vertical (z) velocities as well as the average veloc- 
ity magnitude in the salt provided an indication of the flow stability in the fuel salt, as well as 
when it became steady. The maximum and minimum values were very ‘spiky’ when gravity 
was applied. Cases that diverged during the development of the solution strategy were seen 
to produce sustained unphysically large velocities. 

Heat transfer: The integral of the heat flux across the cladding outer surface and to the inner sur- 


face of the duct provided an assessment of the convergence of the heat transfer solution, 
and the former also allowed energy conservation to be judged, because all of the heat source 
must cross the cladding to the coolant. The expected value is known from the fuel salt vol- 
ume and volumetric source term, and so can be compared to an absolute value, and not just 
a measure of change in the solution. 

Salt mass: The effect of changing from Boussinesq to ‘full’ buoyancy can be seen in changes in 
the mass of fuel salt, giving a positive indication that the change has worked as expected. 
The stabilisation of this value reflects the stabilisation of the mean fuel temperature. 


A selection of these monitors are plotted in Figure 3.8 for the forced circulation conditions. 


The solution was configured to be considered converged and automatically stopped (further itera- 
tions terminated) when the residuals for all fields fell below 10~° (10~’ for energy) and the relative 
change in all of the four monitors underlined above was less than 10~° over the preceding 5 it- 
erations. Iteration numbers in Tables 3.4 and 3.5 marked with an asterisk* stopped automatically 
at these points when their residual and monitor convergence criteria were all met. The solutions 
were enabled to continue to solve for several thousand more iterations, if convergence had been 
not obtained. 


The attention paid to the numerical determination of convergence criteria in these baseline models, 
and to the methods for ensuring solution robustness will be seen later to be key in enabling the 
minimum computation cost and manual intervention to be applied when running a large number of 
automated simulations. 


Interpretation of Results 


The baseline solutions for the forced and natural circulation conditions contain a large quantity 
of data, and despite the apparent good convergence, it is necessary to assess the result for the 
presence of expected features, and for the the absence of artefacts. 


Solid surface fields: The temperature, total heat flux, radiation heat flux and shear stress distribu- 
tions on wall surfaces were assessed, looking for discontinuities, patterns or manifestations 
of the mesh in the results. These patterns could indicate that the resolution is inadequate, or 
the fluid-solid interface is not working correctly, or another solution defect. The y~* field was 
also inspected to confirm that it was less than one. 
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Figure 3.8: Convergence monitors from the forced baseline resolved assembly case so- 
lution. Vertical iteration markers correspond to events in Table 3.4. The minimum vertical 
fuel salt velocity has been inverted to aid comparison. 
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Conservation: The mass and energy conservation was assessed by checking the inlet and outlet 
mass flows and the energy transfers across interfaces, compared to the energy source. 
Cross-section plots: Planes through the fuel, coolant, cladding and duct regions were used to plot 
key quantities as contour plots at a range of heights. These plots allow the results to be 
interpreted for their meaning as a solution result, but also allow solution artefacts related 
to mesh or poor convergence to be seen. It is also possible to check overall features, such 
as whether the solution is symmetric as expected, that the effect of gravity is in the correct 
direction, or that fluid properties or turbulence fields have values in the expected ranges. 

Line plots: The same quantities as investigated with cross-section plots were sampled on lines 
drawn from the centre to the inter-wrapper gap, through the fuel, coolant, cladding and duct, 
at a range of heights. Similar assessments were made as for the contour plots by assessing 
line plots of these quantities — the visualisation as a line plot identifies features that are harder 
to notice on a coloured contour plot, and vice-versa. 


Examples of the cross-section contour plots at the mid-height of the fuel pin for a range of variables 
are shown in Figure 3.9 for the forced circulation baseline case and in Figure 3.10 for the natural 
circulation case. For the forced circulation case, it is not useful to visualise the temperature and 
velocity fields on the same contour plot for fuel and coolant simultaneously. The fuel is too hot and 
slow, and the scale of the colour bar means that both cannot be distinguished. Line plots allow that 
comparison to be made, however, and further assessment of these results using them are shown 
in Section 3.5 in comparison to a single subchannel model. 


The most obvious result visible in both cases is that there is not a significant variation over the fuel 
pins except at the outer row adjacent to the duct, where the velocities and temperatures change 
rapidly. This result can also be seen in the line plots in Section 4.1 when they are compared 
to an equivalent porous simulation. No variation was seen in the fuel flow pattern that indicated 
widespread asymmetry in the natural circulation. Visualising temperature in the coolant for the 
forced circulation conditions is difficult because the thermal boundary layer is thin compared to the 
velocity boundary layer, as a result of the Pr > 1 coolant salt. 


Turbulence Modelling in the Fuel Pin 


For the natural circulation cases, the models become effectively laminar when converged in both 
the coolant and fuel — the turbulence decayed to almost nothing. This confirms that the damping 
terms in the RANS model used are able to represent laminar conditions. The switch to a truly 
laminar solution changes very little in the solution (and hence the model converges after only 51 
more iterations, Table 3.5). 


It was expected that, once convergence had been achieved, the fuel salt flow in the forced circu- 
lation case could also be made laminar by using Fluent’s ‘laminar zone’ option for the fuel tube 
regions only. However, it was found that it was necessary to keep the RANS model active in the 
fuel; running with the fuel as a laminar region produced erratic maximum salt velocities observ- 
able in the monitors, approximately 10 times higher than with a turbulence model active, and high 
solution residuals were seen. 


From a simple pipe flow Reynolds number estimate using the pin dimensions and the expected 
velocities (0.05 m/s), and from the comparison in Section 3.2, it was expected that the fuel circula- 
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Figure 3.9: Cross-sections of the resolved assembly model under forced circulation con- 
ditions at z = 0.8 m (mid-height). 
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Figure 3.10: Cross-sections of the resolved assembly model under natural circulation 
conditions at z = 0.8 m (mid-height). 
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tion would be laminar. Further investigation described in Section 3.4 showed that the need to use 
a turbulence model was not just a numerical solution artefact, however. The following facts indicate 
that the solution was legitimately picking up on the need to model the presence of turbulence in 
the fuel: 


* The natural circulation cases were able to return laminar flow even with the RANS model 
active means that it would not spuriously ‘insist’ on turbulence being present. 


* The poor behaviour when attempting to run with laminar flow in the salt for the forced circu- 
lation case. 


« When solving with RANS turbulence modelling active in the fuel, the turbulent viscosity ratios 
for the forced circulation case are predicted to be >10 (Figure 3.9f) which indicates non- 
negligible turbulent effects. These ratios are actually higher than in the coolant for the same 
case, although the turbulence intensity in the fuel is only ~1%, whereas in the coolant it is 5 
to 10%. 


Single Fuel Pin 


The validation comparisons described in Section 3.2 ran 


Assessment of 
RANS models in 
buoyancy 
affected vertical 
pipe flow 


successfully with laminar settings in the closed tube do- 
main, as did the low heat load natural circulation cases 
above for the resolved fuel assembly, but the high power 
forced circulation cases did not. 


Validation of Sixth sector 


buoyant flow in assembly 
: : fo closed tube with model resolving 
To understand the difference, and determine if it was Heer generate nite cladding: 


fluid fuel and duct 


present for numerical or physical reasons, a simplified 


version of a whole single fuel pin (with the same dimen- 
Full length fuel 
pin to assess 
buoyancy 
induced 
turbulence 


sions as the pins in the fuel assembly) was created. The 
fluid properties were made constant (sampled at 1200 K) 
the Boussinesq approximation was applied to maintain 
constant density, a uniform heat source of 100 MW/m? 
was used and a constant temperature external boundary condition of 900 K was applied on the 
outside of a cladding. This results in a value of S = 4.78 x 10+. Turning the RANS model (k-w 
SST with low Re corrections) on and off in Fluent repeatably and reversibly generated erratic high 
velocities and high residuals when laminar, which disappeared and converged and a steady, sym- 
metric convection pattern appeared with the turbulence model active. Running the same setup in 
CFX resulted in the same behaviour, which suggested that the behaviour was more likely to be the 
result of a physical turbulent process, not a numerical quirk. 


This single fuel pin model was investigated in detail, and it allowed the conclusion to be drawn that 
the flow of fuel salt for the forced circulation case is expected to be turbulent, and is able to be 
represented by the chosen RANS models. The steps used, and the process for understanding this 
buoyancy driven effect, arriving at this conclusion and creating robust supporting justifications are 
discussed in Appendix A. 


The experiments discussed in Section 3.2 were performed in wider tubes of similar length, with 
lower heat sources, but these changes substantially cancelled out to produce a similar value for the 
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non-dimensional parameter S. S is a Rayleigh number multiplied by r/L, the radius to length ratio 
of the tube, which is 1/50 in the experiments, and (dy,e)/2)/L = 1/340.4 for the fuel pin. Therefore, 
for a condition of a given S, the value of Ra in the fuel pin must be higher. The equivalent value of 
S in the experiments by Murgatroyd and Watson (1970) at the Ra used in the single fuel pin is 


(dfyer/2)/L 


Arex iy) 
. ( 1/50 


) = 325 x 10° 

which is in the range where instabilities, and towards the transition to turbulence, were observed in 
the experiment. In strongly buoyancy affected flows, it is often the case that a critical Rayleigh num- 
ber determines the transition to turbulence, and so the higher Ra in the fuel pin is likely generate 
transition at lower values of S. 


This is an example of how important it is to understand the meaning and validity of the parameter 
range for any validation data. While the fuel pin is within the tested range of S values, it is not in 
the experimental r/L, hence Ra range, and there is not a parameter elucidated by Watson (1971) 
that characterises laminar-turbulent transition explicitly. 


Single Coolant Subchannel 


The resolved assembly model demonstrated that the be- 
haviour in the fuel and coolant can be well represented 


2D slice with 
no flow to 
determine 

effective planar 

conductivty 


by a single subchannel and surrounding three fuel pins. 
A mesh was created to represent this subset of the 
cross-section, but with the same axial extent and res- 


Sixth sector 
assembly using 
non-equilibrium 

porous 
medium 


olution as the 1/6" assembly (Figure 3.11). The mesh 
contained 460,000 cells — a factor of approximately 100 
less. Symmetry conditions were applied on the exposed 
sides of the fuel, cladding and coolant. 


Frictional loss and 
cladding to coolant 
heat transfer 
correlations for 
porous model 


Full length 
single coolant 
subchannel, 
with cladding 
and fuel 


The intention of this model is to assess several different 
flow conditions, and then derive friction factor and heat 
transfer correlations from the results. It will also be used 
as the basis for the SA and UQ assessment described in 
Section 4.3. 


The same model setup (with the duct and gap settings removed) and solution strategy as applied to 
the resolved 1/6" assembly model (referred to as the ‘full’ model below for brevity) was able to be 
directly applied to the single subchannel model, and was still found to be necessary to successfully 
achieve a converged solution. 


Comparison to Full Assembly Model 


To assess the agreement between this subchannel model and the full assembly, the same base- 
line forced and natural circulation conditions were applied. The results of this can be seen in Fig- 
ures 3.12 and 3.13 for the forced and natural circulation conditions respectively. The details of the 
fields are easier to see and compare in the more compact image, compared to the contour plots of 
the whole assembly cross-section. There is a horizontal line drawn on Figures 3.12a and 3.13a — 
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Figure 3.11: Single subchannel mesh with surrounding fuel pins. The extended sections 
of coolant at top and bottom (shown) were retained from the 1/6" assembly model. 


Figures 3.14 and 3.15 show the temperature and axial velocity fields along these lines, compared 
to the equivalent results from the centre channel of the full assembly model. The results show the 


key features of the flow: 


The temperature in the fuel varies significantly across the distance from the cladding to the 
centreline of the fuel pin. 


In the fuel, the velocities are upwards in the centre of the pin, and downwards at the edge. 


In the forced circulation case, the fuel temperatures are much higher than in the coolant, and 
the fuel velocities are much lower. In the natural circulation case they are both closer in scale. 


In the forced circulation case, the coolant flows downwards, with steep momentum and ther- 
mal boundary layers. In the natural circulation case, the coolant flow is upwards and shows 
a laminar velocity profile. 

The ‘dips’ in the middle of Figures 3.14 and 3.15 occur where the flow passes through the 
closest gap between fuel pins. 

The absorption and emission fields for participating thermal radiation are closely related to 
the temperature fields, and have volumetric power densities of 60 to 100 MW/m°, compara- 
ble to the nuclear heat source. 


There are small differences between the profiles along these radial lines, but Figures 3.16a to 3.16c 
show that this is because of accumulated differences in the axial profiles. The radial comparison 
is made to the centre of the full assembly, which will have some variation in temperature and flow 


distribution across its cross-section with height. 
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Figure 3.12: Cross-sections of the single channel model under forced circulation condi- 
tions at z = 0.8 m (mid-height). 
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Figure 3.13: Cross-sections of the single channel model under natural circulation condi- 
tions at z = 0.8 m (mid-height). 
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Figure 3.14: Temperature and axial velocity in single subchannel model compared to the 
centre of the full assembly model for baseline forced circulation conditions. 
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Figure 3.15: Temperature and axial velocity in single subchannel model compared to the 
centre of the full assembly model for baseline natural circulation conditions. 
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c: Fuel axial velocity at the fuel pin centreline. 


Figure 3.16: Comparisons of vertical profiles for the full assembly model and equivalents 
in the single channel model for the forced (left) and natural (right) circulation cases. 
Central fuel pin and nearest coolant subchannel used for full assembly. Note the direction 
of coolant flow. 
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* The variation of coolant and fuel temperature in the forced circulation case is almost identical 
(Figure 3.16a). For the natural circulation case, the conditions at the inlet (base) diverge 
slightly with height (the largest discrepancy is at top/outlet). 


¢ The divergence in temperature in the natural circulation case is caused by a difference in 
coolant velocity (Figure 3.16b). The coolant velocity increases with height in the natural cir- 
culation case — this is caused by a subtle concentration of flow towards the full assembly 
centre, driven by buoyancy, and the increasing radial distance from the duct. There is a sim- 
ilar effect in the forced circulation case, but it is proportionally smaller. 


* The steep change in velocity along the centreline of the coolant channel at a height of 1.6m 
(at the top of fuel pins, Figure 3.16b) for the forced circulation case is not a result of devel- 
oping flow (noting that the flow inlet is at 1.7m). It is associated with the sudden change in 
wall heat flux (from zero) at the start of the solid cladding section, which results in a change 
in turbulence, flattening the coolant velocity profile, reducing the centreline value. 


¢ The variation of fuel centreline velocity is driven by the local coolant temperature on the 
other side of the cladding, and so there is bigger difference between the full and single case 
for natural circulation conditions, although it is relatively small (Figure 3.16c). The presence 
of the inflection in the velocity field near the base of the fuel pin in the forced circulation 
case, and its absence for the natural circulation case is consistent with the cases described 
in Section 3.2 at the equivalent values of S. Similarly, the larger fuel axial velocities in the 
forced circulation case compared to the natural circulation case were expected, given the 
higher buoyant driving forces. 


These differences are relatively small, and perfect agreement between the centre of the full case 
and the single cases is not needed to allow the latter to be used to derive friction factor and heat 
transfer correlations over a range of flow conditions. 


Variation in Flow and Derivation of Correlations 


The three forced circulation conditions from CH1 to CH3 and the natural circulation conditions from 
the three instants in CH2 defined in Table 2.2 were used to provide a range of flows and heat trans- 
fer to enable the range of applicability of friction factor and heat transfer correlations to be broader. 
These flow conditions were applied to the single subchannel model. The CH1 and CH3 flows were 
continued from a partially converged baseline (CH2) result, this provides a saving in computation 
cost, but was also found necessary to obtain convergence. The solution strategy developed in Sec- 
tion 3.3.3 is therefore not universal, and would need some modification to successfully start these 
cases with different flow conditions from an initialised state. The same was found with the natural 
circulation conditions at 150s and 500s — they were started from the 250s solution. 


Correlations® for friction factor and Nusselt number in the coolant were required, parameterised 
by Prandtl and Reynolds numbers, so that the dimensionless quantities can link the flow and heat 
transfer processes with the material properties of the salt, giving a more transferable result. The 
characterisation of convection is discussed in Volume 2 (Section 2.1.2). 

Because the temperature and heat flux fields vary along the length of the coolant channel, and there are coolant density and 
viscosity changes associated with this, the flow is never completely ‘fully developed’ with constant heat flux (as idealised in 


mathematical analysis). However the differences are expected to be small, given the relatively slow changes in conditions 
and the small hydraulic diameter relative to the length. 
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Planes were created across the coolant and cladding outer surface in the model at 0.2 m inter- 
vals from 0.2m to 1.4m along the active length (the first and last 0.2m were not considered, to 
reduce inlet and end effects). The CFD solution was sampled on these planes® to return a range 
of quantities at each height: 


Coolant channel mass flow rate, W. 


Coolant channel area, A. 


* Mass flux weighted'® average (bulk) coolant temperature, T.. 


Area weighted coolant density, p. 


Area weighted cladding outer surface (forming the coolant flow path) wall temperature, Ty. 


Area weighted cladding wall shear stress, 7. 


Area weighted cladding wall total heat flux, g, (including the net radiative heat transfer from 
the cladding surface to the coolant). 


¢ The volumetric heat source in the fuel, qvo/. 


¢ The fuel centreline temperature at the given height, 7<;. 


The cladding inner surface, 7;. 


The bulk temperature, 7., was used to evaluate the coolant material properties that are used in 
further calculations. The HTC, h, was obtained from the cladding wall heat flux 


and converted to Nusselt number using 


—hDp 


N 
ee 


The Reynolds number at each height was calculated using the coolant mass flow and subchannel 


hydraulic diameter and area 
_ WD» 


Re= 
‘e Ai 


and the Darcy friction factor was obtained from the wall shear stress 


f= Ar 
~ W2/(2pA?) 


These results are plotted in Figure 3.17 and 3.18 for the forced and natural circulation cases re- 
spectively. 


Results extracted on planes were effective in this case because the flow was steady and the domain was a simple 1D 
channel. For more complex cases it may be preferable to average fluid properties over a short section of volume, and 
average wall properties over the appropriate adjacent bounding surface. Care is needed in selecting the volumes and 
surfaces, and the mass flux weighted average temperature still needs to be calculated, not a density or volume weighted 
average. 

° surface-integrals/mass-weighted-avg in Fluent — using a different surface average will give the wrong result. 
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Figure 3.17: The variation of Nusselt number and friction factor with Reynolds number for 
the three forced circulation conditions in the single coolant subchannel case. A relevant 
value that would be predicted for turbulent flow from the literature is shown, along with 
the best fit equation delivered from the CFD results. 
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Figure 3.18: The variation of Nusselt number and friction factor with Reynolds number for 
the three natural circulation conditions in the single coolant subchannel case. Relevant 
values of Nu and f that would be predicted for laminar flow from the literature is shown, 
along with the best fit values delivered from the CFD results. 
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Forced circulation: For turbulent flow, Nu increases with Re, and also depends on the Pr of 
the fluid. Friction factor decreases with Re at transition Re values, and is constant at high Re. 
Figure 3.17 shows that these trends are observed in the single channel results. The variation in 
Re within each of the three region results (CH1 to CH3) is caused by the variation in viscosity with 
temperature. 


For Nu, Kakag et al. (1987, Chapter 7, for a fluid with Pr = 10 and Re > 10*) suggests that the 
expected value in a triangular rod array with p/d = 1.2 should be approximately 1.1 times greater 
than the Nu in a circular tube. The expected circular tube Nu is represented by the Gnielinski corre- 
lation (Volume 6, Section 3.1, also available in Incropera et a/., 2011 or Rohsenow ef al., 1998). The 
same reference suggests that the expected friction factor for this array would be approximately the 
same as would be found for a circular tube, although this guidance is for Re = 10*. The Haaland 
correlation (Volume 6, Section 3.1, Kakag et al., 1987 or Rohsenow et a/., 1998) for circular pipes 
(applying zero roughness) is included for reference. The Gnielinski correlation requires a friction 
factor; this was provided by the values of f found from the CFD. 


Curve fits have been applied to the CFD results, with the expressions and coefficients shown on 
the figure. The Nu expression has the same form as the Dittus-Boelter correlation (producing the 
same power of = 0.8 for Re, and = 1/3 for Pr as found in several other correlations), and the 
friction factor expression has the same form as the Blasius correlation. 


The results obtained from the CFD solution are not in particularly good agreement with the lit- 
erature expectations — the Nu values follow a similar trend but are approximately 1.6 to 2 times 
the Gnielinski values. The friction factor behaviour is similar, being approximately 1.4 times higher 
than the Haaland values. The CFD results (and the resulting best fit correlations) simultaneously 
combine the effects of the subchannel shape, turbulence, buoyancy, thermal radiation and property 
variation; the contributions of each are hard to separate out without rigorous study, that is often not 
possible in an industrial design or assessment process. There are several reasons why the lack of 
agreement could occur: 


The assessments in Section 3.1 indicated that (for a circular tube at Re = 5300) using the 
low-Re corrections to the k- w SST model were likely to over predict Nu and f, but not by as 
much as seen here. It is also the case that in the natural circulation results (below) the laminar 
Nu is substantially larger than the value predicted by the literature, and that is unaffected by 
the turbulence model. 


The transition Re means that the shape of the channel, and hence the velocity profile (in- 
cluding the narrow gap between pins), has more influence, compared to high Re flow where 
the increased mixing causes the velocity to be more uniform. 


The inclusion of the contribution of thermal radiation to the coolant increases the apparent 
convective HTC. It is appropriate to do this because the target application of the correlations 
are porous models and system codes, neither of which would model the radiation separately 
at the pin array scale. 


While the values of Bo shown in Table 2.4 mean that substantial buoyancy effects are not 
expected, the buoyancy opposed flow is expected to enhance heat transfer to an extent. 


Convection correlations generally assume that fluid properties are constant across a cross- 
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section. The high Pr of the coolant and the strong dependence of its viscosity on temperature 
mean that there is approximately a 100 K drop in temperature near to the cladding surface 
(Figures 3.12a and 3.14). This means that there is a localised decrease in viscosity near the 
wall, increasing the gradients and increasing the heat transfer. This is discussed in Volume 
6, Section 3.1, and that guidance is used to estimate the size of the effect. Evaluating the 
coolant viscosity at 850 K, giving w- (bulk) and at 950 K, giving pw (wall) 


w 


(=) = 1.482491 — 1.044 


it can be seen that approximately 4% of the difference between the results is predicted to 
be attributable to the effect of property variation (although the value of n = 0.11 is only 
valid for Re > 104, and the effect at lower values could be different). It will also be shown in 
Section 4.3 that the coolant viscosity has a substantial influence on the cladding temperature. 


Natural circulation: For laminar flow, constant values of Nu are expected, and the friction factor 


is expected to vary as 
fp fRe 


Re 

where fRe = 64 for a circular tube, for example. Figure 3.18 shows that this behaviour is observed 
inthe CFD results. For Nu, at the lower value of Re for each of the flow conditions (corresponding to 
different times in the SBO transient) there is an increased value. This is at the base of the coolant 
channel, where there is a strong inlet effect due to the short distance from the inlet to the base 
of the fuel pins. The expected value of Nu for a constant heat flux boundary for laminar flow in a 
triangular array of p/d = 1.2 is given by Kakag et al. (1987, Chapter 7) and plotted for comparison. 
The same reference provides a value of fRe = 99.8, which is also plotted for comparison. 


Best fit values have been determined from the CFD results and are also shown. The quality of 
the fRe best fit and agreement between the literature value and the CFD is very good. For Nu 
the CFD predicts a value of Nu that is about 1.3 times higher. Compared to the discussion of 
forced circulation conditions above, there are no turbulence modelling effects, and the velocity 
profile should be similar to the literature value (the friction factor agreement supports this). The 
temperature variation with wall distance (Figures 3.13a and 3.15) is smaller, so the effect of varying 
properties will be less. Some enhancement to heat transfer for this buoyancy assisted laminar flow 
is expected, but the effect is expected to be small (Section 2.2.4). The contribution to cladding 
surface heat flux from thermal radiation is proportionally greater in the natural circulation cases 
(approximately 15 to 20% for the natural circulation case compared to 7 to 9% for the forced 
circulation case), and so this will also increase the apparent convective HTC. 


While some details remain that could be pursued further, it is judged that the correlations derived 
are suitable to use as the basis of a porous model representation, as long as the simulations stay 
within or close to the Reynolds number range used to determine them. 
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Solid Fuel Approximation 


Much of the solution difficulty associated with solving these cases is related to the circulating fuel 
salt. It was postulated that an approximation could be made that the fuel was solid, but with an 
enhanced effective conductivity to replicate the effect of the circulating flow. This would allow the 
coolant heat transfer and cladding temperature to still be evaluated with minimal approximation by 
resolving the coolant flow path, and the energy distribution and centreline temperature in the fuel 
to be predicted by a simplified equivalent solid model. 


From the results extracted at each plane, the relationship between the centreline temperature of 
the fuel and cladding was used with an analytical expression for conduction in a cylinder with a 
uniform heat source (Incropera et a/., 2011) to calculate an equivalent conductivity, if the fuel was 


a solid. 
k = Qvoi(dfuet /2)? 
A(T = Ti) 
This value is approximately 5 times higher than the fuel’s conductivity for the forced circulation 


conditions and 2 times higher for the natural circulation conditions. 


By making the fuel a solid with a constant conductivity of ke, evaluated at mid-height (z = 0.8 m) 
then the coolant temperature and velocity can be matched to the circulating fuel case, as can the 
centreline temperature in the mid-part of the fuel. However, the details of the fuel’s radial temper- 
ature distribution, and the centreline temperatures at either end do not match well (they are not 
reported here). The treatment of the fuel as a solid also does not represent a significant com- 
putational cost saving because the coolant flow still needs to be resolved, requiring a substantial 
number cells. 


The insight gained from this assessment is that there is little net axial heat transfer by the cir- 
culating fuel salt over the majority of the fuel pin height. The local energy source and the local 
cladding surface heat flux were closely balanced — the circulation acts to reduce the fuel centreline 
temperatures but does not redistribute the energy source. This is not true at either end of the fuel 
pin, where the velocities reduce and reverse, and the distribution of energy is more complex. This 
means that any simplified approximation of the internal fuel behaviour details (including a porous 
representation) is likely to contain inaccuracies at the top and bottom of the assembly. This could be 
overcome by creating a spatially varying effective conductivity, but deriving this, and generalising it 
to other flow conditions is a non-trivial task, and may not be able to be easily linked to parameters 
with physical significance. 


An alternative approach to using an auxiliary evaluation to predict the internal state of a fuel rod, 
without solving for the flow in it explicitly, is shown in Section 4.2. 
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3.6 


Performing the Analysis UO 


Transverse Effective Conductivity 


The single subchannel models described above provide 


2D slice with 
no flow to 
determine 

effective planar 

conductivty 


the coolant pressure loss and heat transfer characteris- 
tics from the cladding to the coolant. Beyond this, the 


Radiation contribution 
to planar effective 
conductivity 


remaining properties needed for a porous model depend 
on the geometry and properties of the materials only, ex- 


Sixth sector 
assembly using 
non-equilibrium 

porous 
medium 


cept for the net heat transfer in the planar (transverse, 


or cross-bundle) direction. The effective conductivity in 
these directions in the pin bundle is the net result of the 
combination of several simultaneous heat transfer mech- 


Full length 
single coolant 
subchannel, 
with cladding 
and fuel 


anisms: 


* Conduction through the coolant, fuel and cladding. 
* Thermal radiation transport through the fuel and 
coolant. 


These mechanisms and their inter-relationships are too complex to be evaluated reliably by simpli- 
fied calculations, and so are simulated using a 2D slice model. This allows a temperature depen- 
dent material property to be derived to represent them in a porous model. 


The fuel, coolant and cladding mesh from all pins in the 1/6'" assembly model at a single plane 
were used, bounded by the inner surface of the duct, which was held at a known temperature using 
a large HTC. The salt and cladding were all created as solids, with their normal material properties, 
except the fuel salt, which was given an effective conductivity 5 times normal, as derived above, to 
account for the increased cross-pin heat transfer caused by the circulation"'. 


A heat source was applied to the fuel that was 1000 times lower than the baseline natural circulation 
case. The intention is to assess the shape of the temperature profile from the bundle centre to 
the duct and the difference between centre and duct in as close to isothermal a condition as 
possible. Using this, the fixed duct temperature can be varied to evaluate the dependence of the 
shape and temperature difference (hence effective conductivity) on array temperature, given the 
T* dependence of the thermal radiation contribution. This is an approach that has been used 
previously for modelling spent PWR fuel in dry storage, where thermal radiation transport plays an 
important cooling role (US DOE, 1996). 


Four thermal radiation conditions were considered: 


Transparent: DO radiation, but with zero absorption. 

Low absorption: DO radiation, with an absorption coefficient of k = 300 m~?. 

High absorption: DO radiation, with a higher absorption coefficient of k = 3000m7?. 
No radiation: The high absorption limit, where no radiation transport is modelled. 


An assessment of the optical thicknesses for the high and low absorption case is given in Table 2.3. 
The cases with radiation used a refractive index of 1.3 and a surface emissivity of 0.8, as used in 
the 3D models. The wide range of conditions from transparent to opaque is assessed because 


11 Testing the derivation of planar conductivity without this enhancement showed that it added only a relatively small contribu- 


tion to the result, so the use of a single representative value is an acceptable approximation. 
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c: Various homogeneous conductivities. 


Figure 3.19: Temperature profiles from the centre of the pin bundle slice to the corner of 
the duct (left) and to the middle of the duct wall (right). 
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the absorption coefficient in the salts is unknown and variable. The approximation has been made 
for simplicity here that the fuel and coolant absorption coefficient are the same. The effect of 
independent variation in the absorption coefficients is assessed in Section 4.3. 


Results for the low and high absorption cases are shown in Figure 3.19a and 3.19b for five duct 
temperatures from 800 to 1200 K each. The results are presented relative to duct wall temperature 
Tw for each temperature. Two results are plotted on each figure, one for the temperature profile 
from the centre of the bundle to the corner of the duct (along a symmetry plane) and the other in 
the x-direction from the centre to the middle of the flat face of the duct wall. Using these two results 
provides a cross-check and assurance that the representation of the effective conductivity does 
not exhibit a preferred direction when compared to a homogeneous conductivity (see below). The 
changes in temperature gradient with distance can be associated with passing through the various 
materials, noting that the heat source is only present in the fuel pin'*. There is a clear non-linear de- 
pendence of the centreline-to-duct temperature difference on duct temperature. These differences 
are substantially larger, for a given duct temperature, in the high absorption case compared to the 
lower absorption case, indicating that the effective conductivity is lower with more absorption. 


The relationship between these simulations that contain all of the planar heat transfer mechanisms 
and an effective single value of conductivity at each temperature is established by creating a series 
of cases where the coolant, cladding, and fuel conductivities are all assigned the same constant 
value, and the same total heat source is distributed evenly through all three materials. Conductivity 
values from 0.5W/m K to 4W/m K were chosen and the results are shown in Figure 3.19c. The 
overall shape of the curves agrees with the inhomogeneous models for both the centre-to-corner 
and centre-to-duct results if overlaid, but the structure of the pins is not visible in the results because 
of the uniform conductivity and heat source. The range of centre to wall temperature differences 
covers the range seen in the inhomogeneous models. 


The results for centreline to wall temperature differences vs homogenous conductivity are plotted 
in Figure 3.20 and a best fit is found with the following form 


kerf = a(T _ Tw)? +C 


The best fit coefficients (a, b and c) and the relative fit error for the duct corner and flat results are 
shown in the figure legend. The corner and duct variants are effectively the same. 


This relationship allows an effective conductivity to be evaluated at each duct temperature for the 
inhomogeneous models based on the predicted centreline to duct temperature. These relation- 
ships are plotted in Figure 3.21, where the transparent setup has the highest effective conductivity 
and greatest dependence on temperature, with increasing absorption producing lower and flatter 
curves. A third order polynomial (see Volume 2, Section 3.4.3.3) was fit to each result to allow its 
easy inclusion in Fluent. 


12 The differences in the details of the shapes of the temperature profiles between the x = 300 m~! and 3000 m~? results is 
caused by repeated radiation absorption and (isotropic) re-emission, which allows the radiation to transport energy around 
the fuel pins through the coolant ‘beyond line-of-sight’, as well as by reflection and absorption/emission from the cladding 
surfaces. This is most prominent in the intermediate optical thickness k = 300m! case. This is relevant in the full 
assembly case near to the duct, where there are high temperature gradients over the outer row of pins. 


61 of 103 


Performing the Analysis 


Although these relationships were derived using an approximately isothermal whole bundle cross- 
section, they are also applicable to porous model cases where there are significant localised pin- 
to-pin temperature gradients, caused by large planar heat fluxes — the local effective planar con- 
ductivity would depend on temperature in the vicinity of each pin. 
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Figure 3.20: Relationship and curve fit between centreline to wall temperature difference 
and homogeneous conductivity. 
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Figure 3.21: Effective planar conductivity as a function of duct temperature for the four 
considered radiation conditions. 
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The models described in Section 3 and the information derived from them can now be used for two 
purposes: 


¢ Implementing a porous representation of a fuel assembly that is suitable for inclusion in a 
reactor scale simulation. 


¢ Exploring the effect of uncertainty in material properties on quantities of interest. 


The rest of this section shows example of how to accomplish these applications. However, before 
describing these activities, which would be the point at which this modelling would start to be 
applied to making real decisions, it is worth recapitulating some of the caveats and limitations on 
what has been derived so far. 


¢ Wire wraps have been ignored in this modelling for simplicity. Real fuel assemblies will have 
wire wraps or spacer grids as part of their construction, which will need to be accounted 
for. There will also be inlet and outlet structures (filters, nozzles and fuel supports) that will 
contribute to the pressure loss across the core, that would also need to be present in a whole 
assembly model. 


* The position reached regarding the appropriate choice of RANS turbulence model could 
benefit from further investigation and refinement. For example, the intermittency model 
showed promise, but would need further assessment of its robustness and generality to 
make use of it. 


« Rigorous mesh sensitivity studies have not been undertaken, and would be a necessary part 
of making use of the results to make decisions. Examples of mesh sensitivity studies are 
described in Study A and Study D. 


These points notwithstanding, the purpose of the case study has not been the determination of a 
specific model and its numerical results. Instead, its intention has been to demonstrate the process 
of planning and executing the simulations, extracting the results and using them to derive charac- 
teristics of the fuel bundle, and then the implementation of those characteristics, with comparison 
of their predictions with a reference solution. 
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4.1 


1 


Application of the Results 


Porous Assembly 


To test the implementation of a porous assembly model, 


the same baseline forced and natural circulation con- 
ditions as evaluated by the resolved assembly model 
in Section 3.3 were assessed. By carefully aligning the 
mass flow rates and heat source between the two mod- 
els, a detailed comparison can be made of the quality of 
the agreement between the approaches. 


The porous 1/6" assembly model covers the same re- 
gion of space as the resolved model. As shown in Fig- 
ure 4.1, the domain was split into three segments axially: 


* A 10mm lower region below the rod bundle. 


2D slice with 
no flow to 
determine 

effective planar 

conductivty 


Radiation contribution 
to planar effective 
conductivity 


Sixth sector 
assembly using 
non-equilibrium 
porous 
medium 


Frictional loss and 
cladding to coolant 
heat transfer 
correlations for 
porous model 


Full length 
single coolant 
subchannel, 
with cladding 
and fuel 


* The 1600 mm rod bundle region. 


* A 100mm upper region above the rod bundle. 


and into four regions radially: 


ri: 


r2: 


3: 


gap: 


The majority of the cross-section from the centre to the mid-line of the final row of pins. 
This region contains a heat source in the rod bundle axial region only, and has the same 
geometrical characteristics (porosity, hydraulic diameter and surface area density, influencing 
pressure loss and heat transfer) as the subchannel models. 

From the mid-line of the final row of fuel pins to the outer edge of the final row. This region 
also contains a heat source in the rod bundle axial region, but the change in coolant path ge- 
ometry cause by being adjacent to the duct results in a different porosity, hydraulic diameter 
and surface area density. 

From the outer edge of the final row of fuel pins to the duct. This region contains no heat 
source, and has different geometrical characteristics to r1. The r2 and r3 regions form a ‘wall 
subchannel’ between the centres of adjacent fuel pins and the duct, the characteristics of 
which! are available in the literature (Kakac ef al., 1987, Chapter 3). 

The duct is represented by a thin two sided wall so the gap region represents the combination 
of the duct and coolant in the gap. There is no heat source, and the geometrical properties 
are those of the gap, not the bundle. 


The mesh contains 46,000 cells for the fluid, with the same number of solid cells co-incident with 
them (an identical overlaid mesh) to represent the solid parts of the region?. This is necessary for 


the non-equilibrium porous approach, where the solid can have its own temperature field that is 


different to the flow fluid temperature. This requires the surface area density (interface surface area 
per unit volume), Ax, of the fluid-solid interface and the associated HTC, hg, to be specified. 


The axial splits are required to place the heat source in the correct places, and the gap requires 


a separate region because it is not connected to the bundle. However, splitting the bundle into r1 


to r3 is not essential — it was performed to test various options for modelling the behaviour of the 


Compared to the ‘central subchannel’ assessed in Section 3.5. There is a third ‘corner subchannel’ type at the duct corner 
that has also been analysed in the literature, but the distinction has not been considered here. 
2 92,000 cells in total, compared to 44 million cells in the resolved model, a factor of approximately 480 less. 
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Separate bundle regions: 
r1: covering all of pins to 
the mid-line of the final row. 
r2: from the mid-line of the 
final row to the edge of the 
final row. 

r3: from the edge of the 
final row to the duct. 


Duct and gap region. 


The duct is a thin (zero 
thickness) two sided wall. 


Comparison of the location 
and mesh of resolved pins 
to the porous model. 


Significantly reduced mesh 
resolution in the planar and 
axial directions compared to 
the resolved sixth assembly. 


Mesh regions. 


Figure 4.1: Geometry and mesh of the 1/6" porous fuel assembly. 
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coolant in the near duct region. This followed the descriptions of Wang et al. (2020), who were 


concerned with duct temperatures and inter-wrapper heat transfer in a LMFR application. They 


used a separate zone between the rod bundle and the duct too, that treated the equivalent of r3 as 
a pure fluid layer (e.g. without any porous modelling). 


The details of the chosen setup for the porous model were determined by a series of precursor 


simulations with representative losses and heat transfer. The modelling choices and findings were: 


The physical velocity formulation was assessed as well as the default superficial velocity 
formulation. The guidance in the solver documentation (ANSYS, 2020b), suggests that the 
physical formulation is able to resolve velocity gradients more accurately. The physical ve- 
locity, Up, in a porous region is related to the superficial velocity, U;, by the porosity, ¢, (the 
fraction of the volume that the fluid occupies) by 


U, = Up, 


Regardless of which formulation is used, however, Fluent requires that inlet velocities are su- 
perficial, and that the porous loss coefficients are specified as being relative to the superficial 
velocity. 


Attempting to run with r3 represented by a clear fluid (no solid) or using the physically realistic 
higher porosity of r2 and r3 was problematic when using the physical velocity formulation. The 
flow tended to traverse across the bundle at the inlet and concentrate in r3. The topic of the 
physics and modelling of flow parallel to a porous/non-porous interface is well described in 
the literature (for example Ochoa-Tapia and Whitaker, 1995, de Lemos, 2012 or Al-Aabidy 
et al., 2020) but involve a relatively complex implementation. It is possible that using these 
methods and resolving the mesh in the radial direction in r2 and r3, to properly capture the 
shear and turbulence would be successful, but this would be counterproductive to the aim of 
minimising cell count and model complexity, so they have not been employed here. 


A slip (zero-shear) condition was applied to both sides of the duct wall. The friction associated 
with the duct is included in the porous resistance term already, and there are not enough cells 
to resolve the near-wall gradients. 


Applying a viscosity scaling function (such as the Brinkman correction) was found to make 
little difference, and was omitted. 


From these initial investigations, and the fact that each additional per-region setting adds to the 


complexity of including them at large scale in a reactor model, a relatively simple approach was 


taken. For the forced circulation case: 


The r1 region was represented using the subchannel values of ¢, D;, and A;,, and the f and 
Nu correlations derived from the subchannel model in Section 3.5. 


Two successful configurations for modelling r2 and r3 were found: 
1. Using the physical velocity formulation for all regions and with r2 and r3 represented 
using their own values of Dy, and Ax, but the same value of ¢ as r1. 
2. An alternative configuration, using the superficial velocity formulation for all regions and 
the physically accurate value of ¢ for r2 and r3 (as well as their own values of D, and 
Ats): 
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The same f and Nu correlation as r1 was used in each configuration. Using the correct value 
of ¢ in r2 and r3 with the physical velocity formulation resulted in unphysical discontinuities 
in the turbulence and velocity fields at the r1 to r2 interface. This is largely caused by the 
fact that by default Fluent does not account for the presence of the porous medium in the 
modelling of turbulence at all. An option to avoid this is to run the model with a laminar setup 
(removing all turbulence); this was successful in removing the discontinuities, and the results 
were similar to, but not quite as good as, the ‘alternative’ configuration above, which used 
the superficial velocity formulation. 


* The gap region used its own value of ¢, D,, As, and the same f and Nu correlation as r1 to 
r3 for the turbulent flow conditions. 


For the natural circulation case: 


« The r1 region was represented using the subchannel values of ¢, D, and A,, and the fRe 
and Nu value derived from the subchannel model in Section 3.5. 


¢ The r2 and r3 regions were represented using their own values of ¢, D, and A;,, and the 
same fRe and Nu values as r1. 


* The gap region used its own value of ¢, Dp, As; and fRe = 96 and Nu = 140/17, which are 
analytical results for laminar flow in an infinite parallel plate channel with constant heat flux 
(Kakag et al., 1987, Chapter 3). 


Some of the approximations made are likely to lead to some reduced accuracy of the flow and heat 
transfer in the near-duct region. If the details of the near-duct region were to become of increased 
importance and interest, the modelling details in r2, r3 and the gap could be revisited. 


Porous Model Implementation 


More precursor calculations and modelling parameter evaluations were necessary to create a 
porous model in CFD than explicitly resolving the geometry. These are areas where both con- 
ceptual and mathematical errors are likely, and so require particular attention in verification. 


Geometrical Characteristics: The expressions and values for ¢, D;, and A;, for the four regions 
are given in Table 4.1. The wall subchannel expressions for r2 and r3 use w = (2/3)h, which is the 
distance covered by r2 and r3 (centre of the last pin to the duct, Section 2.2.1). The gap geometry 
has been calculated by treating it as an infinite parallel plate channel, which ignores corner effects. 


Boundary Conditions and Source Term: The inlet and outlet conditions for the porous model 
are the same as the resolved model, except that the velocity specified is the superficial velocity 
across the entire inlet face, giving for example 0.2597 m/s compared to 0.6638 m/s (Table 2.2). 
Care is needed with the flow areas to ensure that exactly the same inlet mass flow rate as the 
resolved model is introduced. 


The heat source was distributed evenly across the solid regions of r1 and r2, using the same axial 
power profile for scaling with height as the resolved case. Care was also needed with the details of 
the cross-sectional areas of the regions to introduce exactly the same power input as the resolved 
case. It should be noted that the source terms are related to the geometrical volume of the region, 
and are not scaled by 1 — ¢ in the solid. 
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r1 r2 and r3 gap 
Ar V3p?/4—1d?/8  (w—d/2)p—nd?/8 
Ao V3p?/4 (w — d/2)p 
Pr nd/2 md/2+ p 
p Ar /Ao As /Ao Weap/(Wegap + 2tauct ) 
Dr 4A; / Pr 4As / Pe 2Wegap 
Afs Pr /Ao Pr /Ao 2/(Wgap + 2tduct) 
¢ 0.37021 0.52766 0.60784 
Dy (m) 0.0058783 0.0063330 0.0062 
Ags (m2/m3) 251.92 333.28 392.16 


Table 4.1: Expressions for the geometrical characteristics of regions and their resulting 
numerical values. A is the cross-sectional area available for coolant flow in a subchan- 
nel and Aj is the overall cross-sectional area of the subchannel, including the region 
occupied by cladding or fuel. 


Equivalent Material Properties: The coolant properties can be used unmodified in the porous 
model and in the gap, as Fluent will correctly weight the contributions to density, specific heat 
capacity and thermal conductivity from the coolant and the HT9 material of the duct by ¢ and 
1 — ¢ respectively. If modelling the heat transfer across the inter-wrapper gap becomes more of a 
focus, then the effective conductivity in the planar direction would need to be modified compared 
to the axial conductivity to represent the layout of the duct and coolant gap (as shown in Volume 2, 
Figure 3.5). 


The solid in the bundle region is intended to represent the combined effect of the fuel salt and 
cladding as an equivalent solid material. The conductivity is split into an axial component and a 
planar component represented using an orthotropic material. 


The effective axial conductivity is evaluated from material properties only, so, as a simplification, 
no contribution from the natural circulation of the fuel is included. 


Ketad Aclad ZG kfuel A fuel 
Aclad a Afuel 


Kequiv axial = 


It should be noted that the k properties here are temperature dependent polynomials, and kag 
is a piecewise function, so the resulting equivalent property must be algebraically or numerically 
derived to be a piecewise polynomial in T too. 


The effective planar conductivity, including thermal radiation, ker, has already been derived in Sec- 
tion 3.6, and the function derived for x = 300 m~! was used to align with the baseline conditions. 
This function already accounts for the fact that the cladding and fuel only occupy 1 — ¢ of the 
cross-section, and contains the contribution from the molecular conductivity of the coolant. How- 
ever, Fluent will perform a porosity weighting of the solid property that would ‘double count’ the fluid 
contribution. This is anticipated by specifying the conductivity of the solid in the planar directions 


to be 


kerf _ Pk coolant 
Kequiv,planar = ~ 1-6 


It should be noted that this expression is strictly only valid for an equilibrium porous model (where 
there is no difference in local fluid and solid temperature) and it is also only ever an approximation to 
segregate the fluid and solid effects; conduction and thermal radiation in the coolant will be altered 
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by the velocity and temperature profiles when there is flow, including by the additional transport 
created by turbulence. 


Calculating effective values for density and specific heat capacity have a further subtlety — it is 
not possible to specify a temperature dependent density for a solid in Fluent (because no thermal 
expansion is modelled), however, it is only the product pc, that matters (and only during a transient 
simulation). The equivalent density is evaluated as 


p 2. Pclad Aclad on Pfuel A fuel 
equiv — 
Aelad ar Afuel 


and a nominal value of this, p,, is established at a representative temperature (say, 1000 K). This 
is used as the material property in Fluent, and also to derive a modified equivalent specific heat 
capacity, which can vary with temperature when implemented 


Cc a PcladCp,clad Aclad a Pfuel Cp, fuel A fuel 
p,equiv — 
PnAclad =P PnAfuel 


This evaluates to the correct temperature varying function as ppCp,equiy because p, is cancelled 
out. This is also needs to be derived as a piecewise polynomial because of the piecewise expres- 
sion for the HT9 specific heat capacity. 


Pressure Loss: The resistance (pressure loss) caused by the porous medium is implemented 
as a momentum source term in Fluent in terms of superficial velocity U.: 


3 3 
$;=- s DjjpUs; + a Ci 5 \Us| Us 


j=l j=l 


Where S; is the momentum source (N/m?) in direction i, Dj; is the viscous loss tensor (m~2), and 
C;; is the inertial loss tensor (m~+). The directions used are the Cartesian axes, with the axial loss in 
the z direction, and the planar loss coefficients in the x and y directions, so there are no off-diagonal 
terms (D and C are diagonal). Only the axial losses are of importance in the fuel assembly (there 
is negligible cross-flow), and the transverse direction losses would not be substantially different, 
so only a single value of either D or C needs to specified for all three directions (making the 
approximation that the losses are isotropic). These assumptions and approximations will not always 
be valid, depending on the geometry and application/flow. 


Comparing this equation to the expression for pressure loss per unit length in a channel flow (also 
with units N/m?) 

aP _ HR) wy, 
This equation is framed in terms of the physical velocity, which is also used to evaluate the Reynolds 
number dependent friction factor, as derived in Section 3.5. For the turbulent forced circulation 
cases, only the inertial term will be used, so an expression for C is required; for the laminar natural 
circulation cases, only the viscous term will be used, so an expression for D is required. More 


complex loss characteristics could use both viscous and inertial terms simultaneously. 


Equating the expressions for the inertial term and relating the superficial and physical velocities 
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using U; = #U, gives 


1 
C=p|U;|Us = 


Algebraically expanding the form of f(Re,) further is not helpful because the expression cannot be 
directly used in Fluent, and to retain the Reynolds number dependence derived in Section 3.5, it 
must be implemented as a User Defined Function (UDF). Similarly equating the viscous term 


f (Rep) 
h 


1 
Dpus, = 5 1Up|Up 


expanding the laminar friction factor expression f = fRe/Re, and expanding Re = p|U,|Dp,/u 


fRe 1 
DwU, = —p|U,|U 
Lb ReD, 2° P| Pp 
fRew F 
2p|Up|Dj, 
fRe 


203 


U, 
Dy = Uelee 


This is a single number, where the viscosity term has been cancelled, so can be implemented in 
the Fluent interface without a UDF. 


Heat Transfer: The heat transfer between the fluid and solid at any point in the porous region is 
a volumetric energy source term 
Qvol = Asshes(T¢ = T;) 


and 


Because k varies with temperature, hs cannot be specified as a single number and requires a 
UDF, even when Nu is assumed constant (in the natural circulation cases). Where Nu depends on 
the local value of Pr and Re (in the forced circulation case), then hs, must be specified as a UDF 
implementing the relationships derived in Section 3.5. 


The solid contains the heat source, which at every location is transferred to the fluid via the HTC 
and surface area density, by a temperature difference. From this equation, it can be seen that 
the interpretation of the solid temperature, 7, in the non-equilibrium porous formulation is that it 
is the cladding outer surface temperature necessary to dissipate the internal source. There is no 
representation of the ‘inside’ the solids. This is only conceptually accurate at steady-state, and 
could potentially be significantly inaccurate in a transient analysis for two reasons: 


¢ The porous solid temperature is assumed to be representative of all of the material in the 
cross-section of a fuel pin with an equivalent heat capacity, pc,. From the results of the re- 
solved rod cases (compared in more detail below) the internals of the fuel are significantly 
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hotter than the surface because they contain the heat source. There is a quantity of energy 
that would need to be accumulated or lost in a transient that represents that elevated tem- 
perature, and this would affect the timescale of the transient response. The presence of this 
‘stored’ energy would not be accounted for using the setup described here. 


* The uniform temperature of the solid implies that it has an infinitely high conductivity (zero 
Biot number) at the scale of individual pins, because there are no local/internal temperature 
gradients represented. In a resolved model, if the coolant flow was suddenly increased, then 
the fuel pin surface would drop rapidly, but the internal temperature profile of the fuel would 
not respond instantaneously. In a porous model, the thermal mass of the whole pin cross- 
section would notionally need to be cooled to reduce its temperature, which would happen 
more slowly. Therefore the transient solid temperature would no-longer represent the local 
surface condition, and the predicted rate of heat transfer to the coolant would be inaccurate. 


Modelling this more accurately within the porous model framework would require a combination 
of modified equivalent material properties for pc, and tracking an additional variable (or variables) 
representing the internal state (temperature profile) of the pins. 


Porous Model Solution: Similar initialisation and numerical settings were applied to the porous 
model as to the resolved model, and some of the same incremental solution strategy steps were 
used, but less caution was needed. It was found that in the natural circulation case, starting the 
solution with turbulence active and only making it laminar after the initial solution stage was a 
reliable approach. Conversely, it was found that when using the superficial velocity formulation (in 
the ‘alternative’ forced circulation configuration), this condition was harder to converge successfully 
and a period of laminar solution at the start of the simulation before activating turbulence modelling 
was successful in this case. 


Comparing Porous and Resolved Assembly Models 


The porous and resolved full assembly model have deliberately been configured to have the same 
mass flow rate and heat source. Therefore, assessing the agreement between the models in 
coolant velocity and temperature does not provide significant insight into the accuracy of the porous 
modelling approach, but it does allow a verification cross-check that the porous implementation is 
correct. Comparisons of temperature prediction of the solid in the porous medium to the cladding 
outer surface in the full model, and an assessment of the overall pressure drop do test the accuracy 
of the approach. 


Correctly predicting overall pressure drop is important because it will determine the flow in each 
fuel assembly in a reactor scale model in response to the above and below core conditions, which 
is likely to be the most important result. The static pressure difference, dP in the coolant from top 
to bottom of the fuel pin section has two contributions: 


dP = dP hydro + AP sriction 


The hydrostatic pressure gradient increases the pressure with height (and so dPhydro is negative). 
Flow from top to bottom in the forced circulation case leads to a positive pressure difference for the 
frictional loss term, dP Fiction. The natural circulation case, with flow from bottom to top, gives a neg- 
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ative dP fiction. The pressure differences and their contributions are shown in Table 4.2, comparing 


the full resolved assembly, single subchannel and porous model. 


Full Single Porous Porous (alt) 
Forced circulation 
dP hydro -45197 -45204 -45176 -45132 
OP friction 12574 12814 13751 12479 
dP -32623 -32390 -31425 -32651 
Natural circulation 
dP hydro -45442 -45440 -45368 
AP friction -490 -468 -590 
dP -45931 -45908 -45958 


Table 4.2: Comparison of pressure difference (N/m*) contributions for different model 
implementations. 


These results are visualised in Figure 4.2a, along with the vertical profile of coolant density and 


temperature, comparing the porous to the full assembly. 


The hydrostatic part of the pressure difference compares well between the full and porous 
models. It is the integral of the density field with height (Figure 4.2b), and this density field 
arises from the coolant temperature field (Figure 4.2c). 


The coolant temperature field for the forced circulation case agrees well over the fuel pin 
height. For the natural circulation case, the gradient with height is similar, but has an offset 
that appears at the base where the coolant enters. The porous model does not represent the 
entrance effects in the coolant flow paths, nor the fuel pin recirculation end effects that occur 
in the full model, leading to a different local heat transfer. 


The natural circulation case has a very low frictional loss — this loss balances the differences 
in hydrostatic driving force produced by temperature differences in the rest of the reactor to 
produce natural circulation. 


The dP for the ‘alternative’ configuration for the forced circulation case (using the superficial 
velocity formulation and accurate value of ¢ in r2 and r3) is also shown in Table 4.2 where 
it can be seen that it produces good agreement with the full assembly. However, it produces 
worse agreement with the flow distribution across the cross-section and with the coolant and 
solid temperatures (these are not reported). 


Cross-sections of the temperature and velocity fields are shown in Figure 4.3, which can be com- 


pared their equivalents in Section 3.3.4, but the porous results produce little visible structure to the 


flow. The per-cell solution values are shown instead of smooth fields (which are interpolated by the 


post-processor) to emphasise the coarseness of the mesh, and also because interpolation of such 


coarse results can introduce misleading visual artefacts. 


Radial line profiles of the velocity and temperature fields produce more meaningful results, and 


results at three heights are shown in Figures 4.4 and 4.5 for the forced and natural circulation 
cases respectively. 
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Figure 4.2: Comparisons of centre subchannel of the full assembly model to the centre 
of the porous model for the forced (left) and natural (right) circulation cases. 
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a: Coolant temperature (K) for forced circulation. b: Vertical velocity (m/s) for forced circulation.. 
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c: Coolant temperature (K) for natural circulation. d: Vertical velocity (m/s) for natural circulation. 


Figure 4.3: Cross-sections of coolant temperature and vertical velocity in the porous 
assembly model at z = 0.8 m (mid-height). 


* The agreement with the coolant velocity is relatively good across the cross-section, and the 
locations near to the duct where it is either increased or reduced compared to the main bun- 
dle in the full model are captured, although a qualitative comparison has not been made. The 
velocities are comparable because the physical velocity formulation is used; for the super- 
ficial formulation, the porous model results would need to be divided by their local value of 
¢. 

* The agreement with coolant temperature in the main part of the bundle and the reduction 
toward the duct is well captured, but the extent of the reduction at the duct is under-estimated 
in the porous fluid. 


« For the natural circulation case, the agreement between the solid temperature in the porous 
model and the full model is good, and the reduction in pin temperature towards the duct is 
represented. The interface between the coolant and the outer surface of the cladding in the 
full model is marked with crosses for each pin. The effect of the combination of conduction 
in the fuel cladding and radiation can be seen in the full model, where there is a small 
temperature offset between the inward facing side of each pin compared to the outward 
facing side. The difference at each height between the porous fluid and solid temperature 
indicates that the choice of the non-equilibrium approach was successful. 
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* For the forced circulation case, the porous solid temperature is generally in good agreement 
with the full model — the differences are larger at higher values of z in the bundle (closer to the 
coolant inlet). Contributions to the transverse thermal transport from turbulence are lower in 
the porous model and likely to introduce some inaccuracy. Turbulence is not correctly affected 
by the presence of the solid, and is under-estimated, resulting in lower thermal diffusivity, 
including in the near duct region. 


* The ‘alternative’ forced circulation configuration (not shown), exhibits lower turbulence in the 
bundle. Its coolant velocity prediction in the main bundle section is lower and does not overlap 
with the full model in each coolant passage, and the solid temperatures are higher. It does, 
however, produce better dP agreement. 


The relatively good comparison of cladding outer temperature is reassuring that the heat transfer 
processes have been characterised correctly, but they do not allow any prediction of the internal 
state of the fuel pin. An approach to obtaining such a prediction is described in Section 4.2 below. 


These comparisons demonstrate the process of ‘closing the loop’, where the resolved model has 
been evaluated, a porous model created, and then compared back to the resolved model as a 
reference. More improvements to the porous model could be made to improve the agreement, 
especially in the duct region if local temperatures are of interest, or to tune the loss representation 
to produce closely matched pressure drops. 


This process could be done in an ad-hoc manner, or could be automated using a tool such as 
Dakota to perform an optimised calibration. However, at each possible iteration of model refine- 
ment, the question should be asked about how much additional effort is justified? How much more 
time and cost in porous model development, with potentially increased complexity in setup of a 
reactor scale model, is worth how much extra accuracy? Particularly when weighed against other 
sources of inaccuracy or uncertainty. 
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Figure 4.4: Comparison of resolved and porous assembly radial profiles of temperature 


and axial velocity at three heights for forced circulation case. 
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Figure 4.5: Comparison of resolved and porous assembly radial profiles of temperature 


and axial velocity at three heights for natural circulation case. 
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Surrogate Model of Fuel Conditions 


The porous model described above gives a good pre- 
diction of cladding temperatures over the whole assem- 


Full length 
single coolant 
subchannel, 
with cladding 
and fuel 


bly, but no internal details of the fuel. It is often the case 
that only the maximum values are of interest. This sec- 
tion demonstrates the construction of a surrogate model 


Automation of 
subchannel for 
SA, UQ and 
surrogate 
models 


where, using the single subchannel case (Section 3.5), 
the flow rate, decay heat and inlet temperature will be 
varied over values relevant to the SBO conditions, anda 
model representing the maximum fuel salt temperature 
will be derived. 


Given the assumption of quasi-steady conditions® then the surrogate model could be used as an 
auxiliary calculation to a porous or system code model to post-process the expected maximum 
temperature at any instant, or it could be used on its own in a design optimisation. 


The Dakota toolkit was used to automate the selection of the cases to be run and also generated 
the surrogate model — a quadratic polynomial response surface model was chosen. 


* A polynomial response surface model is expected to perform well with the smooth variation 
of output quantity, but may give inaccurate results at extremes of the variable range, and is 
liable to be invalid for extrapolation outside of them. 


* A Gaussian Process model could also have been chosen, and would have the advantage of 
providing an estimate of the uncertainty at each point. However, it would be more complex 
to evaluate, and the polynomial model performance is sufficiently accurate to make the per- 
location uncertainty estimate unnecessary. 


More details on the characteristics of surrogate models can be found in Volume 4 (Section 4.4). 
The natural circulation case setup was used, and three input variables were chosen to be varied: 


* The volumetric heat source in the fuel, varying uniformly from 1.8 to 3 MW/m?. 
* The inlet velocity to the coolant, varying uniformly from 0.03 to 0.065 m/s. 
¢ The coolant inlet temperature, varying uniformly from 830 to 860 K. 


The definition of these inputs was provided to Dakota, and 64 evaluations requested, chosen by 
Latin Hypercube Sampling (LHS). It was necessary to write some interface functionality to sup- 
ply these variables to placeholder locations in a Fluent input journal, launch the Fluent run, and 
then extract and return the quantity of interest (fuel salt maximum temperature) to Dakota. Python 
scripting was used for the interface. 


The subchannel model is well suited to this kind of automation because, based on the solution ap- 
proach developed for it in Section 3.5, it already had a suitably parametrised run journal, and solved 
robustly, stopping with a self identified definition of convergence. This is beneficial because it al- 
lows solutions that need more iterations to converge to have them, but also to not waste computing 


i.e. the prediction remains valid as long as none of the inputs change too rapidly, although the definition of what ‘too rapidly’ 
is, is not addressed here. 
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time for those that do not. In addition, it means that all cases are converged to an equivalent state 
of accuracy. The model is also fast to run, taking 20 to 30 minutes to converge on 16 cores. This 
level of automation could equally well be applied to the porous model, but could be prohibitively 
expensive to apply to the resolved assembly case. 


The presence or absence of known ‘sentinels’ (specific phrases in run log file) was used to detect 
run failure and non convergence. Solution robustness is necessary for the combinations of inputs 
provided to Fluent by Dakota, but is hard to ensure a priori across the whole parameter space. 
Therefore, if a solution fails to pass its automated converge checks, it can be manually inspected 
to assess if the solution is still acceptable. In some instances, a CFD solution may diverge or hang, 
and will only run successfully with some changes to the solution strategy, or computer failures may 
interrupt the solution. In these cases, Dakota can be restarted part-way through an interrupted 
assessment, without needing to repeat the already completed CFD evaluations. 


The form of the polynomial model returned by Dakota is a sum of the coefficients, C, each multi- 
plying the product of the input variables, x, raised to various powers, p 


Ny N; 
f(x) = pee Te" 
k=1 i=1 


In the surrogate model derived for this example N, = 10, for the various combinations of powers 
of the N; = 3 input variables. The coefficients and corresponding powers are shown in Table 4.3. 
Dakota also returns the quality metrics for the surrogate model: 


RMS = 0.3388 K 

mean absolute error = 0.2723 K 
max absolute error = 1.151 K 
R? = 0.9995 


This indicates that the fit quality was very good, and that the surrogate can be used, within the 
range of the input variables, to provide a rapid evaluation of what the single coolant subchannel 


CFD model would have predicted, with an error typically less than +1 K. 


Cx Pk. = Pk,2— Pk,3 
2.925205 x 102 0 0 0 
9.152148 10-5 1 0 0 
—1.299638x 10-12 2 0 0 
—3.188494x 10-4 1 1 0 
—5.110499x 10-8 1 0 1 
—2.110934x 103 0 1 0 
1.829253 x 104 0 2 0 
3.868871x107! 0 1 1 
3.630437 «10-4 0 0 1 
3.995353 10-4 0 0 2 


Table 4.3: Coefficients, C, and p (powers of x) for the example maximum fuel tempera- 
ture surrogate model. 


The order of the input variables corresponds to the order that they were supplied to Dakota, so as 
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an example using the 250s (baseline) conditions, for an input vector of 


x, = 2.40812 x 10° W/m? 
X2 = 0.040058 m/s 
x3 = 842.22K 

f(x) = 917.98 K 


The surrogate prediction compares favourably with the CFD result of 917.53 K. 


This process produces the example maximum fuel temperature surrogate model, but it also gen- 
erated 64 CFD solutions, sampling a wide range of the relevant parameter space. Many quantities 
of interest can be extracted and surrogate models fit for them simultaneously, maximising the ben- 
efit derived from the CFD runs performed. The data generated could also be used for a range of 
other purposes, either by Dakota, or by further manual processing. For example, it could be used 
to revisit and improve the range of applicability of the heat transfer and friction factor correlations. 


Uncertainty Quantification of Single Coolant Channel 


The process of arriving at the comparison between the 
porous and resolved models identified potential signifi- 


Full length 
single coolant 
subchannel, 
with cladding 
and fuel 


cant sources of uncertainty: 


+ If simple correlations from the literature had been 
used to create the porous model inputs, then signif- 


Automation of 
subchannel for 
SA, UQ and 
surrogate 
models 


icant differences to the CFD prediction would have 
arisen. 


However, there are uncertainties in the CFD re- 
sults, particularly related to the ability of RANS tur- 
bulence models to represent the fuel and coolant 
flow. 


The porous implementation of the same conditions as the resolved case does not produce 
perfect agreement because of the limitations and simplifications inherent in the porous ap- 
proximation. 


The importance of these is hard to judge until they can be compared to the dependence of the 
results of interest to safety or operational limits, and to the differences caused by other sources of 
uncertainty. There is potentially substantial uncertainty present in the material properties. The final 
application of the models derived in this case study is to use the Dakota toolkit and the solution 
automation mechanisms described above to evaluate the effect of variability in salt properties. 


The type of UQ that Dakota performs is ‘forward propagation of parametric uncertainties’, where in- 
dividual high-fidelity model evaluations are based on sampled distributions of input values, and the 
outputs are assembled using statical methods to provide meaningful insight. The Dakota training 
materials* describe A Practical Process for UQ comprising seven points: 


4 dakota.sandia.gov/training/dakota-training-materials 
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1. Determine your UQ analysis goal: 


« What are the key model responses (quantities of interest)? 
¢ What kinds of statistics or metrics do you want on them? 


ie) 


. Identify potentially influential uncertain input parameters, which includes parameters that 
influence the trend in a response, as well as those that influence variability in a response. 


3. Characterise input uncertainties and map them into Dakota variable specifications. 


4. What are the model characteristics/behaviours? Simulation cost, model robustness, input/output 


properties such as kinks, discontinuities, multi-modal, noise, disparate regimes. 
5. Select a method appropriate to the variables, goal, and problem. 
6. Set up Dakota input file and interface to simulation. 


7. Run study and interpret the results. 


The goals for the MSR fuel assembly UQ (Point 1) are to determine the influence of uncertainty 
and variability in salt material properties on three quantities of interest: 


Max clad T: the maximum temperature in the fuel pin cladding solid®. 
Max fuel 7: the maximum temperature in the fuel salt. 
Max U,: the maximum axial velocity in the circulating fuel salt. 


The probability distribution of the range of predicted values that are possible for these quantities 
is sought, as well as the contribution to variability in each from the chosen uncertain properties. 
This will provide the contribution of these variables to the expected uncertainty band that should be 
applied to predictions of these quantities of interest. It will also allow a judgement to be made about 
which properties should be subject to further modelling or measurement investigation. These three 
quantities were already used as solution convergence monitors, so were readily available from 
each CFD solution, and their expected behaviour was already well understood. 


The modelling performed in the proceeding sections addresses Point 4 — the cost and robustness 
of the model is well understood, and kinks and discontinuities are not expected, except at laminar- 
turbulent transition. For this reason, the UQ study is restricted to turbulent flow conditions, and 
the forced circulation case setup for the single coolant subchannel model with baseline (CH2) flow 
conditions was used. 


Points 2 and 3 are discussed in Section 4.3.1 below, Points 5 and 6 are described in Section 4.3.2, 
and Point 7 in Section 4.3.3. 


Selection of Uncertain Input Parameters 


The properties of the intended fuel and coolant salts contain significant uncertainty, and are as- 
sumed to represent a substantial contribution to the overall uncertainty in any assessment. This 
assumption is tested here by performing what would be best considered to be a screening study, 
where relatively wide uncertainties are assumed and both SA and UQ are performed at the same 
time on the same set of CFD simulations. 


5 Abbreviating ‘cladding’ to ‘clad’ for compactness in figures and tables. 
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Not all material properties have been chosen to be varied, and there is no simultaneous varying of 
flow conditions and material properties. Studying a wide range of variables at once is possible using 
the methods described here, but some caution and thought is warranted, rather than attempting a 
brute force approach and varying everything. This is because the complexity of the responses is 
unknown at this initial stage, and the number of runs necessary, as well as the task of assessing 
the adequacy of the outcome and methods, can easily grow to be become difficult to manage and 
interpret. 


The chosen properties, the form of their Probability Density Function (PDF) (normally or uniformly 
distributed) and the assumed uncertainty in them is specified in Table 4.4. The uncertainty in the 
viscosity and thermal conductivity was chosen as the largest values for these properties from any of 
the salts described by Romatoski and Hu (2017). That reference assumes that these uncertainties 
are 20 values (+2 standard deviations, 95% coverage), but acknowledges that there is ambiguity 


about whether the source data is a lo or a 20 value. Because the salts used in this case study are 
less well known, with some of their properties estimated, then the stated uncertainties represent 


+1o, representing one standard deviation (68% coverage) as a ‘standard uncertainty’ (JCGM, 
2008). This means that, for example, if 15% represents 1c, the uncertainty in the value is twice as 
much as if 15% was assumed to represent 2c, i.e. that the PDF is twice the width. 


Normal standard uncertainty 
fuel Op/OT +15% 
fuel k +15 % 
fuel py +20% 
coolant k +15% 
coolant w +20% 
Uniform min max 
fuel « (m~') 100 1000 
coolant « (m~') 100 1000 


Table 4.4: Chosen uncertain variables and their estimated uncertainty. 


When varying a property it was not considered adequate to vary it as a uniform value applied 
everywhere in the fluid in question. The local variation in temperature near to surfaces, and the 
local variation in 4 and k that this produces is expected to lead to differences in heat transfer. This 
is a reason that these properties were selected as a subset of properties, in preference to specific 
heat capacity, because there is no guidance available on the variation of c, with temperature. A 
given perturbation to the » and k properties was implemented by multiplying the entire polynomial 
by a factor and specifying these new polynomial coefficients in Fluent. This altered the way that 
the interface to Dakota was implemented. Dakota specifies a property to be varied as a single 
number, and therefore, for each normally distributed variable, Dakota was requested to select from 
a standard normal distribution (mean = 0, variance = 1), and the corresponding uncertainty value 
and polynomial manipulation was implemented in the python interface, which was developed as a 
‘driver’ for Fluent. 


The error associated with fitting the fifth order polynomial to the exponential viscosity function is 
not accounted for. For the nominal (zero uncertainty) functions the error in this polynomial fit is 
0.43% for the coolant and 0.89% for the fuel (< 1% error being the criterion used to select the 
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4.3.2 


Application of the Results 


necessary polynomial order) and while this is small, it is an example of a range of small sources 
of uncertainties that can be introduced as part of the practical implementation of the modelling 
process that are often overlooked, or not well captured and dealt with rigorously. Completeness 
and comprehensiveness in the treatment of all sources of uncertainty in a complex model is likely 
to be either impossible, or disproportionate when considered as part of a graded approach (Volume 
1, Section 2.2.5). 


In addition to % and k in both coolant and fuel, uncertainty in the temperature coefficient of the 
density of the fuel salt was included. Contradictory evidence exists from some historical references 
(such as that reported in Janz et al., 1975), and so the value of 0p/OT (-1.015kg/m?K as the 
baseline value in Section 2.2.3.3) was varied. This is a different type of polynomial manipulation to 
that applied to ~ and k because instead of shifting a whole curve to be greater or smaller value, 
this changed the slope only. The value of fuel density at 1200 K was used as the reference point, 
with the gradient ‘pivoting’ around this. 


Dakota allows a mixture of forms of PDF for uncertain variables to be assessed, and the normal 
distributions for thermophysical properties were able to be combined with uniform distributions for 
absorption coefficients. The range of absorption coefficients was chosen to span an optically thick 
to an optically thin situation (Table 2.3), and the independent effects of the fuel and coolant salt 
absorption can be established. 


The potential variation in material properties is assumed to have no correlation between them, so 
they can be treated as independent uncertain variables. 


Sensitivity Analysis and Uncertainty Quantification Method 


The chosen approach for SA and UQ of the uncertainty in material properties is to apply a Polyno- 
mial Chaos Expansion (PCE) surrogate model to a sampled set of CFD solutions and to evaluate 
the desired statistical quantities from it. Once the PCE surrogate model has been created, its prop- 
erties allow the statistical moments of the output distribution (mean, variance, skew and kurtosis), 
the linear sensitivity coefficients and Sobol indices to be evaluated analytically. 


Two different methods for creating the PCE were evaluated, each being a cross-check of the other. 


LHS: 128 samples were drawn from the PDFs of the 7 uncertain input variables and a CFD 
model evaluation used for each. The baseline single subchannel converged solution was used 
as a starting point for each evaluation to reduce run-time and reduce the risk of the solution not 
converging in the early phases of establishing the flow. The LHS method is the most common 
sampling method and typically used as the starting point for an SA or UQ study. The samples and 
their solved evaluations can be used for a number of purposes after they have been generated. 
More samples can be added to an existing set if more data is found to be beneficial, but to retain the 
‘latin’ property of covering the parameter space with one sample per ‘bin’, the number of samples 
must double each time. 


The results from the LHS samples themselves can be used directly without applying the PCE, as 
will be shown in Section 4.3.3, and this kind of sampling combined with a straightforward evalu- 
ation of the results would be the only option if some of the uncertain variables were not smooth. 
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Application of the Results 


For example, an LHS sampling approach could vary the RANS turbulence models as one of its 
parameters, but PCE could not, because its inputs must be smooth. 


When building the PCE from the LHS samples, Dakota can perform its own assessment to find the 
best polynomial order for each output variable polynomial, and third to fifth order responses were 
generally selected. The advantage of building the PCE is mainly the additional quantities that can 
be derived from it; to obtain the Sobol indices by LHS sampling alone with 7 input variables would 
require 7 + 2 = 9 times as many CFD evaluations. 


Sparse Grid: The PCE model applied to the LHS samples is a regression fit to that data. An 
alternative method directly uses the orthogonality in relationships between the inputs by creating a 
grid of samples and building polynomials along the various directions. To use the full grid of points 
would need NV? samples for NV uncertain variables with a polynomial order of p. This quickly results 
in a large number of evaluations being necessary. 


An alternative is to apply a ‘sparse grid’ approach, where only a subset of the grid is evaluated. 
Dakota selects the number of evaluations needed based on the ‘level’ of occupancy of the grid. 
For the 7 uncertain variables used here, 137 CFD evaluations were created for sparse grid Level 
2. Level 3 would need 889 evaluations. Similarly to the LHS PCE, the sparse grid output contained 
up to 5th order polynomial responses. 


Both types of PCE were sampled 10° times after being built to provide smooth Cumulative Distri- 
bution Functions (CDFs) with well resolved ‘tails’, in addition to their analytically derived quantities. 
Once built the PCEs could also be potentially re-used as fast running surrogate models for the 
quantities of interest in their own right. 


Uncertainty Quantification Results 


The LHS results were assessed first by using a scatter plot of each of the three quantities of interest 
against the seven uncertain inputs, and by drawing histograms of the outputs — these are shown 
together in Figure 4.6. 


It is good practice to always visualise the results like this to look for features and anomalies in 
the results, rather than relying purely on statistical measures. In particular, the presence of fine 
structure in the response of a variable is significant because it would make the application of a 
method like PCE less likely to be successful. 


For each scatter plot, the partial correlation coefficients and partial rank correlation coefficients are 
shown for the quantity of interest vs. input variable (the rank coefficient is the number shown un- 
derneath). The normally distributed variables are plotted against the value of o chosen by Dakota 
(noting that these are then scaled by the uncertainties in Table 4.4), showing that a span of ap- 


proximately +2o0 was sampled for each. 
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Figure 4.6: Scatter plots and histograms for all three response variables against the seven uncertain variables. 
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Figure 4.7: Example of sparse grid sampling. 


The difference between the random layout of the LHS samples and the sparse grid sampling can 
be seen in Figure 4.7. All 137 samples are plotted on Figure 4.7a (this is a projection onto 2 
dimensions of the 7 dimensional sample space) so the variation of the other parameters occurs 
along lines orthogonal to this plane. 


The standard normal distribution values chosen by Dakota for these variables mainly lie along 
the Oo axis for the other variables (a higher sparse grid level than 2 would place more samples 
towards the corners, to capture more coupled effects, but would require more evaluations), and the 


maximum range is 4.20. 


This is a wider range of variation than the LHS samples, so it was necessary to check that these 
greater variations in properties still produced valid CFD solutions. An example of the response of 
the CFD evaluations on the sparse grid parameter space is shown in Figure 4.7b, where the maxi- 
mum cladding temperature variation is shown against the variation in coolant thermal conductivity. 
The clear trend of reducing cladding temperature with increasing conductivity can be seen, which 
agrees with the correlation coefficients with values that are nearly -1 in the corresponding scatter 
plot. 


The CDFs of the LHS samples are plotted for the resulting variation in the three quantities of 
interest as a result of the variation in all inputs in Figure 4.8. The CDFs from both methods of PCE 
are also plotted, and show effectively identical results as each other, and as the LHS results. The 
advantage of the PCE curves is that they do not contain the same ‘noise’ (variation in gradient from 
point to point) as the LHS, and they can also extend into the low probability high and low response 
‘tails’ of the distribution. 


The mean and standard distribution of the three output responses (which are not normally dis- 
tributed, and have no requirement to be to apply these measures) from the LHS samples corre- 
sponding to the CDF figures are shown in Table 4.5, including the upper and lower 95 % Confidence 
Interval (Cl) estimated by Dakota for each. The mean and o from the sparse grid PCE implemen- 
tation is also shown in this table — the nominal values agree well with the LHS values and lie 
comfortably within the Cls. 
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Max clad T (K) Max fuel T (K) 


Max U, (m/s) 


LHS samples 
lower 95% Cl of mean 990.01 1257.05 0.04060 
mean 991.95 1260.61 0.04106 
upper 95% Cl of mean 993.90 1264.18 0.04151 
lower 95% Cl of o 9.91 18.14 0.00230 
o 11.13 20.37 0.00259 
upper 95% Cl of o 12.67 23.22 0.00295 
Sparse grid PCE 
mean 992.06 1260.61 0.04103 
o 11.39 20.62 0.00257 
Table 4.5: Statistics for UQ for forced circulation conditions. 
Max clad T Max fuel T Max U, 
fuel Op /OT -5.81 2.375 x 10-3 
fuel k -13.8 
fuel uw 9.45 4.061 x 10-4 
coolant k -7.16 -4.95 9.800 x 10-5 
coolant w 7.89 5.56 —1.122 x 10-4 
fuel K 0.0143 3.170 x 10-° 
coolant « -0.00173 


Table 4.6: Linear sensitivity coefficients for SA for forced circulation conditions (omitting small values). 


Max clad T Max fuel T Max U, 

fuel 0p /0T 0.086 0.878 
fuel k 0.492 

fuel w 0.234 0.057 

coolant k 0.471 0.067 0.002 

coolant uw 0.512 0.078 0.002 

fuel K 0.040 0.056 

coolant « 0.005 


Table 4.7: Sobol indices for SA for forced circulation conditions, corresponding to Figure 4.9 (omitting < 10%). 


87 of 103 


Application of the Results 


0.85 4 
£ 
To} 
oO 
6 06> 1 
ou 
(0) 
2 
& 0.45 7 
_) 
& 
) 
0.2; LHS J 
PCE from LHS 
PCE from sparse grid 
0 | 1 1 1 1 


Max clad T (K) 


960 970 980 990 1000 1010 1020 1030 1040 1050 


oO 
(op) 
T 


Cumulative probability 
(eo) 
'N 


LHS 
PCE from LHS 
PCE from sparse grid 


0 
1180 1200 1220 1240 1260 1280 
Max fuel T (K) 


1300 1320 1340 


° ° 
(op) (ec) 
T T 


Cumulative probability 
ic 
'N 


0.2 


LHS 
PCE from LHS 
PCE from sparse grid 


0 
0.032 0.034 0.036 0.038 0.04 0.042 0.044 0.046 0.048 0.05 


Max ws (m/s) 


Figure 4.8: CDFs for the quantities of interest. 
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Application of the Results 


Comparing several methods of calculating the same quantities is recommended, and similar to the 
process of building confidence in a CFD model. Applying SA and UQ successfully also involves a 
gradual accumulation of evidence, insight and decisions, adding understanding of the methods and 
system incrementally — at the early stages of an assessment there is rarely one complete method 
or setup that can be applied without checking its performance and questioning its appropriateness. 


This means that SA and UQ do not ‘come for free’ (even with convenient tools like Dakota). They 
add an additional burden to asking ‘is my physical simulation adequate?’ of also asking: 


* Is my UQ adequate? 

* Have | chosen suitable uncertain variables and distributions? 
* Are the surrogate model polynomials high enough order? 

* Have enough samples been taken? 

* Did the best fits converge well enough? etc. 


If SA and UQ are to be used and relied on to make decisions or justifications, then they need to be 
of high enough quality and sufficiently trustworthy to justify that reliance. 


Sensitivity Analysis Results 


The CDF results discussed above focus on the distribution of the output quantities of interest. The 
same PCE can be used to look at the significance of the uncertainty in each input (distinguishing 
SA from UQ, as discussed in Volume 4). Two results are shown here that allow this: 


Linear sensitivities: The gradient of the response of each quantity of interest to changes in each 

input at the centre of its uncertainty distribution is shown in Table 4.6. Small values have 
been omitted for clarity. These values are from the sparse grid PCE, although the values 
determined for the PCE fit to the LHS samples were very similar. For the normally distributed 
variables, these values are sensitivity per a, for the uniformly distributed radiation absorption 
coefficients, they are per unit k. 
As an example of their application, Figure 4.7b shows that for an increase in coolant k from 
Oo to +1.730 = 26% there is a 10.25 K reduction in cladding temperature in the CFD model 
result. Using the sensitivity coefficient predicts —7.16K/o x 1.730 = —12.4K. The agree- 
ment is reasonable, with the difference arising from the obvious curvature (non-linearity) in 
the actual response. This shows that in this case, projecting the first-order gradient is only 
accurate for small variations. 


Sobol indices: While the sensitivity coefficients can quantify the expected change in an output 
with an input, they are less intuitive to interpret and compare how important the uncertainty 
in each inputs is. The Sobol indices for each input sum to 1, and show the relative contribution 
to variation in the output of each input (although they do not give the sign information that 
sensitivity coefficients do). The Sobol indices® are shown in Table 4.7 (again omitting small 
values) and are plotted in Figure 4.9. 


The Sobol indices give a clear result for each quantity of interest. 


These are the ‘main effect’ indices (see Volume 4, Section 4.5); Dakota also produces the ‘total effect’ indices, which have 
similar values in this case. 
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Max clad 7: Only the coolant thermophysical properties have any significant influence. Variation 
in the fuel properties have no effect, and the coolant absorption coefficient had almost no 
influence. The dependence of the cladding temperature on coolant & was expected, but the 
equal significance of coolant yz was not; this is interpreted as supporting the point discussed 
in Section 3.5.2, where, for this low Re turbulent flow, the effect of the variation of viscosity in 
the near-wall momentum boundary layer is important for heat transfer. 


Max fuel 7: The maximum fuel temperature is affected mainly by the fuel k, and to a lesser extent 
by its viscosity. The influence on the cladding maximum temperature of the coolant pz and k is 
also seen, because this is passed on as an offset to the fuel temperature. The fuel radiation 
absorption coefficient is seen to have only a small effect. 


Max U,: The maximum fuel circulating velocity is sensitive to 0p /OT; this is expected because this 
is effectively the same as 6 that appears in S (Section 2.2.4). Of the other variables, only a 
small sensitivity to fuel % and absorption coefficient was observed. 
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Figure 4.9: Sobol indices for SA for forced circulation conditions, corresponding to Table 4.7. 


The conclusions that can be drawn from the Sobol indices are supported by comparing the rel- 
ative magnitude of the sensitivity coefficients, and also by the correlation coefficients shown on 
Figure 4.6. Both types of correlation coefficient generally agree with each other, and where there 
is a large sensitivity contribution, there is a correlation that is close to +1. Where there is less 


sensitivity, there is less correlation. 


The lack of sensitivity to absorption coefficient is a significant result, because knowing the value 
of this property through all stages of reactor life would be challenging, and modelling its effect in 
detail can be computationally costly and complex. 


The comparative lack of sensitivity of the fuel maximum temperature to 0p/OT is notable, given its 
influence on the maximum fuel velocity. This is consistent with the conclusion reached earlier that, 
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at the mid-length of the pin (where the maximum temperature occurs in the forced circulation case, 
because of the power profile), the heat transfer is mainly in the radial direction in the pin, and so the 
fuel conductivity is the dominant effect, as shown. It is expected that the profile of fuel temperature 
at either end of the pin would be more sensitive to 0p/OT. The interaction of this parameter with the 
turbulence modelling in the fuel salt may also be an area worth investigating further, if desirable. 


This is an initial screening assessment, and it should be noted that although an uncertain input 
may not exhibit sensitivity here, it cannot be excluded from all future consideration based on this 
assessment alone. A different flow condition (such as natural circulation) or a different quantity of 
interest (such as duct temperature or coolant pressure drop) may be significantly more sensitive to 
an input than found here. To be confident in a conclusion, the process of assessing the importance 
of uncertain variables, and the conditions that they are applied in, should be sufficiently extensive, 
and will likely need to be iterative. 
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6 Nomenclature 


Latin Symbols 
A Area, m? 
At Atwood number (At = (p1 — p2)/(e1 + p2)) 
Bi Biot number (Bi = hL/kg) 
Cp, Cy Specific heat at constant pressure or volume, J kg~! K~* 
dorD_ Diameter (D, = 4A;;/Pcs for hydraulic diameter), m 
f Darcy-Weisbach friction factor 
Fo Fourier number (Fo = a,t/L?) 
Gr Grashof number (Gr = gL3Ap/v?p = gL?BAT/v?, using the Boussinesq approxi- 
mation Ap/p ~ —BAT, where AT is often taken as Ty, — Ts,00) 
g Acceleration due to gravity, ms~? 
h_ Specific enthalpy, J kg~+, Heat Transfer Coefficient (HTC), Wm~? K—! or height, m 
I Radiative intensity, W m~? sr~! or Wm~? sr~+ wm? for a spectral density, where sr 
(steradian) is solid angle 
J Radiosity, W m~? 
k Thermal conductivity, Wm-! K-? 
L Length or wall thickness, m 
M_ Molar mass of a species, kg kmol~? 
Ma _ Mach number (Ma = U/a, where a is the speed of sound) 
n_ Refractive index 
Nu Nusselt Number (Nu = AL/k¢) 
p Perimeter, m 
P Pressure (P; = static pressure, Py = total pressure), Nm~? or Pa 
Pe  Péclet number (Pe = RePr) 
Pr Prandtl number (Pr = cpu/ ke) 
q_ Heat flux (rate of heat transfer per unit area, g = Q/A), Wm~? 
Q_ Rate of heat transfer, W 
r Radius, m 
R_ Gas constant (for a particular gas, R = R/M), Jkg~!K~! 
R Universal gas constant, 8314.5 J kmol~+ K~+ 
Rip Thermal resistance, K W~! 
Ra Rayleigh number (Ra = GrPr) 
Re Reynolds number (Re = pUL/,, or for an internal flow Re = WD,,/A.<h) 
Ri Richardson number (Ri = Gr/Re’) 
Sr Strouhal number (Sr = fL/U, where f is frequency) 
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Nomenclature Se, 


«x j<< 


Stanton number (St = Nu/RePr) 

Time, s 

Temperature (7, = static temperature, 77 = total temperature), K 
Wall friction velocity (u, = \/Tw/p), ms~+ 

Velocity, ms~? or thermal transmittance, Wm? K~! 
Specific volume, m* kg~? 

Volume, m? 

Mass flow rate, kg s~+ 
Wall distance, m 


Non-dimensional wall distance (y* = yu, /v) 


Greek Symbols 


Qn SS FE YM KF Hm R WR 


4 


oa 


Thermal diffusivity (a = k/pcp), m?s~? 


Volumetric thermal expansion coefficient (8 = —(1/)(0p/0T)), K~* 
Ratio of specific heats (y = cp/c,) 

Emissivity or surface roughness height, m 

Absorption coefficient, m~? 

Wavelength, m 


Viscosity, kg m~1s~? 


Kinematic viscosity and momentum diffusivity (v = w/p), m?s~? 
Density, kg m~? 

Stefan Boltzmann constant, 5.67 x 10-8 W m-? K~* 

Shear stress, Nm~? 


Porosity or void fraction 


Subscripts and Modifications 


b 
cs 
f 


i 


wn 


8 


~ 


Bulk (mass-averaged) quantity 
Cross-sectional quantity 

Quantity relating to a fluid 

Quantity relating to a particular species 
Total (stagnation) quantity 

Turbulent quantity 

Static quantity or quantity relating to a solid 
Quantity relating to a wall or surface 
Quantity far from a wall or in free-stream 
Average quantity 

Molar quantity 

Varying quantity 
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7 Abbreviations 


CAD 
CDF 
CFD 
CHT 
Cl 
DNP 
DNS 
DO 
EHRS 
FA 
HTC 
LES 
LHS 
LMFR 
MSR 
PCE 
PDF 
PWR 
RANS 
RMS 
SA 
SBO 
SSR-W 
UDF 
UQ 
URF 


Computer Aided Design 
Cumulative Distribution Function 
Computational Fluid Dynamics 
Conjugate Heat Transfer 
Confidence Interval 

Delayed Neutron Precursor 

Direct Numerical Simulation 
Discrete Ordinates 

Emergency Heat Removal System 
Fuel Assembly 

Heat Transfer Coefficient 

Large Eddy Simulation 

Latin Hypercube Sampling 

Liquid Metal-cooled Fast Reactor 
Molten Salt Reactor 

Polynomial Chaos Expansion 
Probability Density Function 
Pressurised Water Reactor 
Reynolds-Averaged Navier-Stokes 
Root Mean Square 

Sensitivity Analysis 

Station Blackout 

Stable Salt Reactor — Wasteburner 
User Defined Function 
Uncertainty Quantification 
Under-Relaxation Factor 
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Single Fuel Pin Turbulence Assessment 


In the 1/6'" fuel assembly model, the fuel pins have high aspect ratio 20 mm long cells for most 
of their length. The single fuel pin mesh used the same resolution as a starting point, but these 
cells were suspected of being problematic, so an additional revised mesh was created, limiting 
the maximum axial cell length to 1mm, so that the cells were not high aspect ratio. This mesh 
contained 3 million cells for the single pin instead of 375,000. 


With this axially resolved mesh, using the simplified conditions described in Section 3.4, it was 
found that when running laminar with a steady-state solver, the maximum velocities were much 
reduced, but the solution residuals would still not converge. 
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Figure A.1: Streamlines for the transient laminar axially resolved single pin case (left) 
where the streamlines are released from the whole cross-section at mid-height and indi- 
cate chaotic flow. For the equivalent case with a RANS model active (right), the stream- 
lines released along a line across the profile show steady, symmetric flow upwards in 
the middle of the pin and downwards at the edge. There are small localised regions of 
significantly faster flow in the transient laminar case than in the turbulent case. 
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Single Fuel Pin Turbulence Assessment 


Starting from a converged RANS steady-state solution, the solution was solved transiently with a 
2 ms time step (keeping Courant numbers at approximately 0.5), and the RANS model was turned 
off. The symmetric flow pattern accelerated, keeping its shape for a period, moving several times 
faster than the initial condition, and then suddenly the maximum velocities dropped and became 
unsteady. The solution residuals were converged at every time step throughout. The interpretation 
is that the large buoyancy forces created the high velocities, generating an instability that eventu- 
ally underwent a process that appeared like laminar-turbulent transition, resulting in chaotic flow 
patterns, that were well resolved by the axially refined mesh (Figure A.1). 


The transient laminar case was solved for 20s and the velocity field was averaged (also gathering 
the RMS values) over the last 11s, once the flow had appeared to reach a statistically stable 
condition. Figures A.2 to A.5 show the centreline and cross-pin profiles of vertical velocity and 
temperature for the transient laminar solution’ and several cases with RANS models: 


k-w SST: the k-w SST model with low Re corrections, run on the coarse (20 mm) axial resolution 
mesh. 

k-w SST high res: the same turbulence model run on the fine (1 mm) axial resolution mesh. The 
results are almost indistinguishable from those for the coarse mesh. A transient version of 
this fine mesh with the RANS model active was solved, and no instabilities occurred — the 
result was identical to the steady-state RANS solutions. 

k-w SST gamma: the ¥ intermittency option was added and the low Re corrections turned off in 
the k-w SST model. This model was started from a converged k-w SST solution and used 
the coarse mesh 

realizable k-«: the realizable k - ¢ model, using the coarse mesh. 


The RANS model results are in reasonable agreement with the averaged transient laminar results: 


* The realizable k - € result for centreline velocity looks convincing in Figure A.2, but Figure A.4 
shows that agreement at the centre is not matched by equally good agreement over the 
cross-section; the shape of the velocity profile does not agree as well with the averaged 
laminar transient result. The velocity and temperature behaviour of this model at the top of 
the centreline also contains features that disagree with the laminar solution. 


* The centreline velocity for the k-w SST based models is higher than the laminar solution 
(thereby predicting a higher fuel salt circulation rate), but the shape of the velocity profile is 
in better agreement than for the realizable k - model. 


* The agreement for temperature is worst for realizable k - e, although not substantially different 
to the k-w SST results — they predict higher temperatures than the laminar solution. The 
exception to this is the intermittency model, which produces a much closer centreline and 
cross-pin profile temperature prediction. 


* The intermittency model prediction for centreline velocity at the base of the pin produces a 
high velocity peak that is not seen in the other models or the laminar results. 


None of the RANS models tested are ideal, and which to choose depends on what other features 
of the analysis of the fuel performance depend on the behaviour. The effect of the difference in 


1 Averaging the flow for 11s was considered sufficient to provide a comparison to the RANS models; a significantly longer 
simulation would be necessary to remove the remaining noise to produce a smooth profile. 
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Single Fuel Pin Turbulence Assessment 


predictions is not limited to the impact on maximum fuel temperatures for thermal reasons. The 
temperatures also impact on neutronic thermal feedback, and hence the reactivity and heat source. 
The circulation rate of the fuel salt in the pin will affect the transport and distribution of DNPs, and 
so could have some impact on the fission power profile. 


While it works well for the temperature prediction in this case, the intermittency model has not been 
selected because it was not designed for application to these conditions, and has a sensitivity to 
solution trajectory as described in Section 3.1. It is possible that it will work well at some locations 
or under some conditions, but be highly inaccurate at others — further work would be necessary to 
prove its robustness. The k-w SST model with low Re corrections was retained in preference to 
the realizable k - ¢ model, primarily because of the features seen in the latter at the top of the tube, 
the slightly better temperature predictions of k-w SST, and the confidence in its ability to revert to 
fully laminar flow where necessary. 


Although the laminar case is a spatially resolved, transient flow with no turbulence modelling ap- 
plied, not enough rigour has been applied to refer to it as a DNS solution. There has been no 
assessment of mesh and time resolution sensitivity, and the result is for a single, simplified case 
without varying density or viscosity. However, the result does offer reassurance that RANS models 
are activating and generating turbulence fields for good reasons, and justifies that it is necessary 
to use one in the resolved fuel assembly for physical, not purely numerical reasons. 


It is surprising to an extent that the RANS models performed as well as they have, because this sit- 
uation is far from the typical flows used to calibrate their parameters. The turbulence generation in 
the models is likely to be partly caused by temperature gradients in the buoyancy production term, 
as well as by shear. The details of the production and low Reynolds number damping processes 
are expected to be important to the behaviour and success of the models. It would be necessary to 
study more features of the result (such as Reynolds stresses) in detail at more locations to ensure 
that the apparent success of the RANS models is for the ‘right’ reasons, and to potentially improve 
or adjust the models for this situation. 


This case also demonstrates a feature of the meshing approach that could cause solution difficul- 
ties in other situations — the very elongated cells in the coarse mesh cannot represent the unstable 
laminar flow. As an instability forms that turns the flow, it cannot be resolved and represented by 
the high aspect ratio cells, so the solution does not perform well. 
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Figure A.2: Single pin centreline vertical velocity. Left: transient laminar case with mean 
value over 11s compared to an instantaneous snapshot, and the envelope of the RMS 
value either side. Right: mean transient laminar solution compared to steady-state RANS 
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Figure A.3: Single pin centreline temperature. Left: transient laminar case with mean 
value over 11s compared to an instantaneous snapshot, and the envelope of the RMS 
value either side. Right: mean transient laminar solution compared to steady-state RANS 
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Figure A.4: Single pin cross-pin vertical velocity profile at mid-height. Left: transient 
laminar case with mean value over 11s compared to an instantaneous snapshot, and the 
envelope of the RMS value either side. Right: mean transient laminar solution compared 
to steady-state RANS solutions. 
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Figure A.5: Single pin cross-pin temperature profile at mid-height. Left: transient laminar 
case with mean value over 11s compared to an instantaneous snapshot, and the enve- 
lope of the RMS value either side. Right: mean transient laminar solution compared to 
steady-state RANS solutions. 
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